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PRINCIPAL NOTATION 


SYMBOLS 

Unless specified otherwise, all variables are nondimensional. 
Symbol Definition 


Cpy Cy 
Cp 

E t 

gc 

hr 

k 

Lr 

M 

n 

Nu N 2 , N z 


P 

Pn 

Q 

R 

R 

Re r 

T 

Uy Vy W 

x>y } z 
y 

S?>, if 


£/ 

k 2 , k 4 
M 

Z, *1 ,C 


Specific heats at constant pressure and volume. 

Static pressure coefficient. 

Total energy per unit volume. 

Proportionality constant in Newton's second law. 

Stagnation enthalpy per unit mass. 

Effective thermal conductivity coefficient. 

Dimensional reference length. 

Mach number. 

Time level. 

Number of grid points in the f , rj, and C directions. 

Static pressure. 

Laminar Prandtl number. 

Vector of dependent variables in the Cartesian coordinate form of the governing 
equations. 

Residual. 

Gas constant. 

Reference Reynolds number. 

Static temperature. 

Velocities in the Cartesian x, y , and z directions. 

Cartesian coordinates. 

Ratio of specific heats, c^/c,. 

Second- and fourth-order explicit artificial viscosity coefficients in constant coeffi- 
cient model. 

Implicit artificial viscosity coefficient. 

Parameters determining type of time differencing used. 

Constants in nonlinear coefficient artificial viscosity model. 

Viscosity coefficient. 

Computational coordinate directions. 

Static density. 

Computational time. 
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SUBSCRIPTS 


Subscript 

Ujt k 
n 
r 
T 

SUPERSCRIPTS 

Superscript 

n 


4 Principal Notation 


Definition 

Denotes grid location in rj, and £ directions. 
Denotes dimensional normalizing condition. 
Denotes dimensional reference condition. 
Denotes total, or stagnation, value. 


Definition 

Denotes time level. 

Overbar; denotes dimensional value. 
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SUMMARY 


A computer code called Proteus has been developed to solve the three-dimensional, Reynolds-averaged, 
unsteady compressible Navier-Stokes equations in strong conservation law form. The objective in this ef- 
fort has been to develop a code for aerospace propulsion applications that is easy to use and easy to modify. 
Code readability, modularity, and documentation have been emphasized. 

The governing equations are written in Cartesian coordinates and transformed into generalized 
nonorthogonal body-fitted coordinates. They are solved by marching in time using a fully-coupled 
alternating-direction-implicit solution procedure with generalized first- or second-order time differencing. 
The boundary conditions are also treated implicitly, and may be steady or unsteady. Spatially periodic 
boundary conditions are also available. All terms, including the diffusion terms, are linearized using 
second-order Taylor series expansions. Turbulence is modeled using either an algebraic or two-equation 
eddy viscosity model. 

The program contains many operating options. The thin-layer or Euler equations may be solved as 
subsets of the Navier-Stokes equations. The energy equation may be eliminated by the assumption of 
constant total enthalpy. Explicit and implicit artificial viscosity may be used to damp pre- and post-shock 
oscillations in supersonic flow and to minimize odd-even decoupling caused by central spatial differencing 
of the convective terms in high Reynolds number flow. Several time step options are available for conver- 
gence acceleration, including a locally variable time step and global time step cycling. Simple Cartesian or 
cylindrical grids may be generated internally by the program. More complex geometries require an ex- 
ternally generated computational coordinate system. 

The documentation is divided into three volumes. Volume 1 is the Anal sis Description, and presents 
the equations and solution procedure used in Proteus. It describes in detail the governing equations, the 
turbulence model, the linearization of the equations and boundary conditions, the time and space differ- 
encing formulas, the ADI solution procedure, and the artificial viscosity models. Volume 2, the current 
volume, is the User's Guide, and contains information needed to run the program. It describes the pro- 
gram's general features, the input and output, the procedure for setting up initial conditions, the computer 
resource requirements, the diagnostic messages that may be generated, the job control language used to run 
the program, and several test cases. Volume 3 is the Programmer's Reference, and contains detailed infor- 
mation useful when modifying the program. It describes the program structure, the Fortran variables stored 
in common blocks, and the details of each subprogram. 

A two-dimensional/axisymmetric version of Proteus code also exists, and was originally released in late 
1989. 
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1.0 INTRODUCTION 


Much of the effort in applied computational fluid dynamics consists of modifying an existing program 
for whatever geometries and flow regimes are of current interest to the researcher. Unfortunately, nearly 
all of the available non-proprietary programs were started as research projects with the emphasis on dem- 
onstrating the numerical algorithm rather than ease of use or ease of modification. The developers usually 
intend to clean up and formally document the program, but the immediate need to extend it to new ge- 
ometries and flow regimes takes precedence. 

The result is often a haphazard collection of poorly written code without any consistent structure. An 
extensively modified program may not even perform as expected under certain combinations of operating 
options. Each new user must invest considerable time and effort in attempting to understand the underlying 
structure of the program if intending to do anything more than run standard test cases with it. The user's 
subsequent modifications further obscure the program structure and therefore make it even more difficult 
for others to understand. 

The Proteus three-dimensional Navier-Stokes computer program is a user-oriented and easily-modifiable 
flow analysis program for aerospace propulsion applications. Readability, modularity, and documentation 
were primary objectives during its development. The entire program was specified, designed, and imple- 
mented in a controlled, systematic manner. Strict programming standards were enforced by immediate peer 
review of code modules; Kemighan and Plauger (1978) provided many useful ideas about consistent pro- 
gramming style. Every subroutine contains an extensive comment section describing the purpose, input 
variables, output variables, and calling sequence of the subroutine. With just three clearly-defined ex- 
ceptions, the entire program is written in ANSI standard Fortran 77 to enhance portability. A master ver- 
sion of the program is maintained and periodically updated with corrections, as well as extensions of general 
interest (e.g., turbulence models.) 

The Proteus program solves the unsteady, compressible, Reynolds-averaged Navier-Stokes equations in 
strong conservation law form. The governing equations are written in Cartesian coordinates and trans- 
formed into generalized nonorthogonal body-fitted coordinates. They are sol\ i by marching in time using 
a fully-coupled altemating-direction-implicit (ADI) scheme with generalized time and space differencing 
(Briley and McDonald, 1977; Beam and Warming, 1978). Turbulence is modeled using either the Baldwin 
and Lomax (1978) algebraic eddy- viscosity model or the Chien (1982) two-equation model. All terms, in- 
cluding the diffusion terms, are linearized using second-order Taylor series expansions. The boundary 
conditions are treated implicitly, and may be steady or unsteady. Spatially periodic boundary conditions 
are also available. 

The program contains many operating options. The thin-layer or Euler equations may be solved as 
subsets of the Navier-Stokes equations. The energy equation may be eliminated by the assumption of 
constant total enthalpy. Explicit and implicit artificial viscosity may be used to damp pre- and post-shock 
oscillations in supersonic flow and to minimize odd-even decoupling caused by central spatial differencing 
of the convective terms in high Reynolds number flow. Several time step options are available for conver- 
gence acceleration, including a locally variable time step and global time step cycling. Simple grids may be 
generated internally by the program; more complex geometries require external grid generation, such as that 
developed by Chen and Schwab (1988). 

The documentation is divided into three volumes. Volume 1 is the Analysis Description, and presents 
the equations and solution procedure used in Proteus . It describes in detail the governing equations, the 
turbulence model, the linearization of the equations and boundary conditions, the time and space differ- 
encing formulas, the ADI solution procedure, and the artificial viscosity models. Volume 2, the current 
volume, is the User's Guide, and contains information needed to run the program. It describes the pro- 
gram's general features, the input and output, the procedure for setting up initial conditions, the computer 
resource requirements, the diagnostic messages that may be generated, the job control language used to run 
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the program, and several test cases. Volume 3 is the Programmer's Reference, and contains detailed infor- 
mation useful when modifying the program. It describes the program structure, the Fortran variables stored 
in common blocks, and the details of each subprogram. 

A two-dimensional/axisymmetric version of Proteus code also exists, and was originally released in late 
1989 (Towne, Schwab, Benson, and Suresh, 1990). 

The authors would like to acknowledge the significant contributions made by their co-workers. Tom 
Benson provided part of the original impetus forlhe development of Proteus, and did the original coding 
of the block tri-diagonal inversion routines. Simon Chen did the original coding of the Baldwin-Lomax 
turbulence model, and consulted in the implementation of the nonlinear coefficient artificial viscosity model. 
William Kunik developed the original code for computing the metrics of the generalized nonorthogonal grid 
transformation. Frank Molls has created a separate diagonalized version of the code. Ambady Suresh did 
the original coding for the second-order time differencing and for the nonlinear coefficient artificial viscosity 
model. These people, along with Dick Cavicchi, Julie Conley, Jason Solbeck, and Pat Zeman, have also 
run many debugging and verification cases. 
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2.0 GENERAL DESCRIPTION 


In this section the basic characteristics and capabilities of the Proteus code are described in general. 
More detailed descriptions can be found in other sections of this manual or in Volumes 1 and 3. 

2.1 ANALYSIS 

Proteus 3-D solves the three-dimensional unsteady compressible Navier-Stokes equations. The 
equations are solved in fully conservative form. As subsets of these equations, options are available to solve 
the Euler equations or the thin-layer Navier-Stokes equations. An option is also available to eliminate the 
energy equation by assuming constant total enthalpy. The governing equations are described in detail in 
Section 2.0 of Volume 1. 

The equations are solved by marching in time using the generalized time differencing of Beam and 
Warming (1978). The method may be either first- or second-order accurate in time, depending on the 
choice of time differencing parameters. Second-order central differencing is used for all spatial derivatives. 
The time and space differencing formulas are presented in Sections 3.0 and 5.0 of Volume 1. Nonlinear 
terms are linearized using second-order Taylor series expansions in time, as described in Section 4.0 of 
Volume 1. The resulting difference equations are solved using an alternating-direction implicit (ADI) 
technique, with Douglas-Gunn type splitting as written by Briley and McDonald (1977). The boundary 
conditions are also treated implicitly. 

Artificial viscosity, or smoothing, is normally added to the solution algorithm to damp pre- and post- 
shock oscillations in supersonic flow, and to prevent odd-even decoupling due to the use of central differ- 
ences in convection-dominated regions of the flow. Implicit smoothing and two types of explicit smoothing 
are available in Proteus . The implicit smoothing is second order with constant coefficients. For the explicit 
smoothing the user may choose a constant coefficient second- and/or fourth-order model (Steger, 1978), 
or a nonlinear coefficient mixed second- and fourth-order model (Jameson, Schrmdt, and Turkel, 1981). 
The nonlinear coefficient model was designed specifically for flow with shock waves. The artificial viscosity 
models are described in detail in Section 8.0 of Volume 1. 

The equations are fully coupled, leading to a system of equations with a block tridiagonal coefficient 
matrix that can be solved using the block matrix version of the Thomas algorithm. Because this algorithm 
is recursive, the source code cannot be vectorized in the ADI sweep direction. However, it is vectorized in 
one of the non-sweep directions, leading to an efficient implementation of the algorithm. The solution al- 
gorithm is described in detail in Section 7.0 of Volume 1. 

2.2 GEOMETRY AND GRID SYSTEM 

The equations solved in Proteus were originally written in a Cartesian coordinate system, then trans- 
formed into a general nonorthogonal computational coordinate system as described in Section 2.3 of Vol- 
ume 1. The code is therefore not limited to any particular type of geometry or coordinate system. The only 
requirement is that body-fitted coordinates must be used. In general, the computational coordinate system 
for a particular geometry must be created by a separate coordinate generation code and stored in an unfor- 
matted file that Proteus can read. However, simple Cartesian and cylindrical coordinate systems are built 
in. 


The equations are solved at grid points that form a computational mesh within this computational co- 
ordinate system. Note that a distinction is being made between the terms computational coordinate system 
and computational mesh . The computational coordinate system refers to the (£,>j,0 system in wdiich the 
governing equations are written. It is determined by supplying a series of points whose Cartesian (xy?,z) 
coordinates are specified, either by reading them from a file or through one of the analytically defined co- 
ordinate systems built into subroutine GEOM. The computational mesh consists of grid points distributed 
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along lines in the computational coordinate directions. These points may differ in number and location 
from those used to determine the computational coordinate system. The number of gnd points in each 
direction in the computational mesh is specified by the user. The location of these grid points can be varied 
by packing them at either or both boundaries in any coordinate direction. The transformation metrics and 
Jacobian are computed using finite differences in a manner consistent with the differencing of the governing 
equations. 

2.3 FLOW AND REFERENCE CONDITIONS 

As stated earlier, the equations solved by Protein are for compressible flow. Incompressible conditions 
can be simulated by running at a Mach number of around 0.1. Lower Mach numbers may lead to nu- 
merical problems. The flow can be laminar or turbulent. The gas constant R is specified by the user, with 
the value for air as the default. The specific heats c„ and c, the molecular viscosity p, and the thermal 
conductivity k can be treated as constants or as functions of temperature. The empirical formulas used to 
relate these properties to temperature are contained in subroutine FTEMP, and can easily be modified if 
necessary. The perfect gas equation of state is used to relate pressure, density, and temperature. This 
equation is contained in subroutine EQSTAT, which could also be easily modified if necessary. All 
equations and variables in the program are nondimensionalized by normalizing values derived from refer- 
ence conditions specified by the user, with values for sea level air as the default. 

2.4 BOUNDARY CONDITIONS 


The easiest way to specify boundary conditions in Proteus is by specifying the type of boundary (i.e., 
no-slip adiabatic wall, subsonic inflow, periodic, etc.). The program will then select an appropriate set of 
conditions for that boundary. For most applications this method should be sufficient. If necessary, how- 
ever, the user may instead set the individual boundary conditions on any or all of the six computational 
boundaries. 

A variety of individual boundary conditions are built into the Proteus code, including: (1) specified val- 
ues and/or gradients of Cartesian velocities u, v, and w, normal velocity K„, coordinate direction velocities 
V(, V v and Vr, pressure p, temperature T, and density p; (2) specified values of total pressure pr, total 
temperature TV, and flow angles; and (3) linear extrapolation. Another useful boundary' condition is a "no 
change from initial condition" option for u, v, w, p, T, p, Pr- and/or 7V- Provision is also made for user- 
written boundary conditions using subroutines BCF and BCFLIN. Specified gradient boundary conditions 
may be in the direction of the coordinate line intersecting the boundary or normal to the boundary, and 
may be computed using two-point or three-point difference formulas. For all of these conditions, the same 
type and value may be applied over the entire boundary' surface, or a point-by-point distribution may be 
specified. Unsteady and time-periodic boundary conditions are allowed when applied over the entire 
boundary. The boundary conditions available in Proteus are described in detail in Section 3.1.7. 

2.5 INITIAL CONDITIONS 


Initial conditions are required throughout the flow field to start the time marching procedure. For un- 
steady flows they should represent a real flow field. A converged steady-state solution from a previous run 
would be a good choice. For steady flows, the ideal initial conditions would represent a real flow field that 
is close to the expected final solution. 

The initial conditions are set up in subroutine INIT. A default version of INIT is built into Proteus that 
specifies uniform flow with constant flow properties everywhere in the flow field. These conditions, of 
course, do represent a solution to the governing equations, and for many problems may help minimize 
starting transients in the time marching procedure. However, realistic initial conditions that are closer to 
the expected final solution should lead to quicker convergence. 

The best choice for initial conditions, therefore, will vary from problem to problem. For this reason 
Proteus does not include a more generalized routine for setting up initial conditions. Instead, the user 
should supply his or her own version of subroutine INIT to set up the initial conditions for the particular 
flow being computed. Details on the Fortran variables to be specified by INIT may be found in Section 
5.1. 
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2.6 TIME STEP SELECTION 


Several different options are available for choosing the time step At, and for modifying it as the solution 
proceeds. At may be specified directly, or through a value of the Courant-Friedrichs-Lewy (CFL) number. 
When specifying a CFL number, the time step At may be either global (i.e., constant in space) based on 
the minimum CFL limit, or local (i.e., varying in space) based on the local CFL limit. For unsteady 
time-accurate flows global values should be used, but for steady flows using local values may lead to faster 
convergence. Options are available to increase or decrease At as the solution proceeds based on the change 
in the dependent variables. An option is also available to cycle At between two values in a logarithmic 
progression over a specified number of time steps. The various time step options are described in detail in 
Section 3.1.9. 

2.7 CONVERGENCE 

Five options are currently available for determining convergence. The user specifies a convergence cri- 
terion e for each of the governing equations. Then, depending on the option chosen, convergence is based 
on: (1) the absolute value of the maximum change in the conservation variables AQ m „ over a single time 
step; (2) the absolute value of the maximum change AQ„„ averaged over a specified number of time steps; 
(3) the L 2 norm of the residual for each equation; (4) the average residual for each equation; or (5) the 
maximum residual for each equation. These criteria are defined in Section 4.1.6. 

2.8 INPUT/OUTPUT 

Input to Proteus is through a series of namelists 1 and, in general, an unformatted file containing the 
computational coordinate system. All of the input parameters have default values and only need to be 
specified by the user if a different value is desired. Reference conditions may be specified in either English 
or SI units. The namelist and coordinate system input are described in Section 3.0. A restart option is also 
available, in which the computational mesh and the initial flow field are read from unformatted restart files 
created during an earlier run. The use of the restart option is described in Sections 3.1.3 and 5.3. 

The standard printed output available in Proteus includes an echo of the input, boundary conditions, 
normalizing and reference conditions, the computed flow field, and convergence information. The user 
controls exactly which flow field parameters are printed, and at which time levels and grid points. Several 
debug options are also available for detailed printout in various parts of the program. The printed output 
is described in Section 4. 1 . 

In addition to the printed output, several unformatted files can be written for various purposes. The first 
is an auxiliary file used for post-processing, usually called a plot file, that can be written at convergence or 
after the last time step if the solution does not converge. Plot files can be written for the NASA Lewis 
plotting program CONTOUR or the NASA Ames plotting program PLOT3D. If PLOT3D is to be used, 
two unformatted files are created, an XYZ file containing the computational mesh and a Q file containing 
the computed flow field. The plot files are described in detail in Section 4.2. Another unformatted file 
written by Proteus contains detailed convergence information. This file is automatically incremented each 
time the solution is checked for convergence, and is used to generate the convergence history printout and 
with Lewis-developed post -processing plotting routines. The contents of the convergence history file are 
presented in Section 4.3. And finally, two unformatted files may be written at the end of a calculation that 
may be used to restart the calculation in a later ran. One of these contains the computational mesh, and 
the other the computed flow field. The contents of the restart files are described in detail in Section 4.4. 

2.9 TURBULENCE MODELS 


For turbulent flow, Proteus solves the Reynolds time-averaged Navier-Stokes equations, with turbulence 
modeled using either the Baldwin and Lomax (1978) algebraic eddy-viscosity model or the Chien (1982) 
two-equation model. Both models are described in detail in Section 9.0 of Volume 1. 


1 It should be noted that namelist input is not part of ANSI standard Fortran 77, but is nevertheless available with 
most Fortran compilers. See Section 2.3.1 of Volume 3 for a discussion of possible computer-dependent features 
in the Proteus code. 
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2.9.1 Baldwin-Lomax Model 


The Baldwin-Lomax model may be applied to either wall-bounded flows or to free turbulent flows. 
For wall-bounded flows, the model is a two-layer model. In the inner region, in addition to the Baldwin- 
Lomax formulation, an alternate expression first presented by Spalding (1961), and later by Kleinstein 
(1967), may be used. For flows in which more than one boundary is a solid surface, the nearest boundary 
is used. The turbulent thermal conductivity coefficient k, is computed using Reynolds analog)'. 

2.9.2 Chien k-t Model 

With the Chien two-equation model, partial differential equations are solved for the turbulent kinetic 
energy k and the turbulent dissipation rate t. These equations are lagged in time and solved separately from 
the Navier-Stokes equations. An LU factorization algorithm is used to solve the k-t equations. 

Since the Chien two-equation model is a low Reynolds number formulation, the k-t equations are 
solved in the near-wall region. No additional approximations are needed. Boundary conditions that may 
be used include: (1) no change from initial or restart conditions for k and e; (2) specified values and/or 
gradients of k and s; and (3) linear extrapolation. Specified gradient boundary conditions are in the direction 
of the coordinate line intersecting the boundary, and may be computed using two-point or three-point dif- 
ference formulas. For all of these conditions, the same type and value may be applied over the entire 
boundary surface, or a point-by-point distribution may be specified. Spatially periodic boundary conditions 
for k and t may also be used. Unsteady boundary conditions are not available for the k-t equations. 
However, unsteady flows can still be computed with the k-t model using the unsteady boundary condition 
option for the mean flow quantities and appropriate boundary conditions for k and t, such as specified 
gradients or linear extrapolation. 

Initial conditions for k and t are required throughout the flow field to start the time marching procedure. 
The best choice for initial conditions for the k-t equations will vary from problem to problem, and the user 
may supply a subroutine, called KEINIT, that sets up the initial values for k and t for the time marching 
procedure. Details on the Fortran variables to be specified by KEINIT may be found in Section 5.1. A 
version of KEINIT is built into Proteus that computes the initial k and t values from a mean initial or re- 
start flow field based on the assumption of local equilibrium (i.e., production equals dissipation.) Variations 
of that scheme have been found to be useful in computing the k-t initial conditions for a variety of turbulent 
flows. 

The time step used in the solution of the k-t equations is normally the same as the time step used for 
the mean flow equations. However, the user can alter the time step for the k-t equations, making it larger 
or smaller than the time step for the mean flow equations, by specifying a multiplication factor. The user 
can also specify the number of k-t iterations per mean flow iteration. 
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3.0 INPUT DESCRIPTION 


The standard input to the three-dimensional version of Proteus consists of a title line and several 
namelists. Additional input may be provided in the form of a pre-stored unformatted file containing the 
computational coordinate system. The calculation can also be started by reading the computational mesh 
and the initial flow field from restart files written during a previous run. This section describes only the 
standard input and the coordinate system file. The restart file contents and format are described in Section 
4.4. 

3.1 STANDARD INPUT 

All of the standard input parameters have default values and do not need to be specified by the user 
unless some other value is desired. The type (real or integer) of the input parameters follows standard 
Fortran convention, unless stated otherwise (i.e., those starting with I, J, K, L, M, or N are integer, and the 
remainder are real.) 2 Note that in most, if not all, implementations of Fortran, namelist names and input 
start in character position 2 or higher in the input line. All of the input, except for namelist IC, is read in 
subroutine INPUT. Namelist IC is read in subroutine INIT. 

3.1.1 Reference and Normalizing Conditions 

Unless specified otherwise, all of the input parameters are specified in nondimensipnal form, with the 
appropriate reference condition as the nondimensionalizing factor. A few words explaining what we mean 
by reference conditions and normalizing conditions, and the differences between them, may be helpful at this 
point. 

The normalizing conditions are, by definition, the conditions used in nondimensionalizing the governing 
equations, and are denoted by an n subscript. (See Section 2.0 of Volume 1.) These normalizing conditions 
are defined by six basic reference conditions, for length, velocity, temperature, density, viscosity, and thermal 
conductivity, which are specified by the user. Reference conditions are denoted by an r subscript. The 
normalizing conditions used in Proteus are listed in Table 3-1. 

Note that for some variables, like pressure, the normalizing condition is dictated by the form of the 
governing equations once the six basic reference conditions are chosen. Unfortunately, some of these may 
not be physically meaningful or convenient for use in setting up input conditions. Therefore, some addi- 
tional reference conditions are defined from the six user-supplied ones. The reference conditions are listed 
in Table 3-2. 

To summarize, the normalizing conditions are used to nondimensionalize the governing equations. The 
average user need not be too concerned about these. The reference conditions are the ones used for 
nondimensionalization of all user- specified input and output parameters. 3 

3.1.2 Title 

TITLE A descriptive title, used on the printed output and in the CONTOUR plot file, up 

to 72 characters long. This is a character variable. 


2 Throughout this User's Guide, elements of the Fortran language, such as input variables and subprogram names, 
are printed in the text using uppercase letters. However, in most implementations, Fortran is case-insensitive. The 
Proteus source code itself is written in lowercase. 

3 Internal to the Proteus computer code itself, variables are generally nondimensionalized by the normalizing con- 
ditions. The reference conditions are used for input and output because they are usually more physically mean- 
ingful for the user. 
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3.1,3 Namelist RSTRT 


The parameters in this namelist control the use of the restart option. The contents of the restart files 
are described in Section 4.4. 

IREST 0 if no restart files are to be read or written. 

1 to write restart files at the end of the calculation. 

2 to read restart files for initial values, and to write restart files at the end of the 
calculation. For turbulent flow, the current calculation must use the same 
turbulence model as the previous calculation. Note that only the initial values 
and the computational mesh are read from restart files. The usual namelist 
input must still be read in. Of course, some input parameters, such as the 
reference conditions or those specifying the grid, must not be changed during 
a restart. 

3 same as the IREST = 2 option, but for cases in which the previous calculation 
was either laminar or used an algebraic turbulence model, and the current 
calculation uses a two-equation turbulence model. 

Two restart files are written and/or read. One contains the computational mesh, 
and the other contains the mean flow field and the k and t values (if a two -equation 
turbulence model is being used.) For IREST — 0 or 1, the initial mean flow field 
will be generated in subroutine I NIT, and the initial k and t values (if necessary) 
will be generated in subroutine KEINIT. The default value is 0. 

NRQIN Unit number for reading the restart flow field. The default value is 1 1. 

NRQOUT Unit number for writing the restart flow field. The default value is 12. 

NRXIN Unit number for reading the restart computational mesh. The default value is 13. 

NRXOUT Unit number for writing the restart computational mesh. The default value is 14. 

3.1.4 Namelist IQ 
n r intout Controls 

The following parameters specify which variables are to be printed, and at what locations in both time 
and space. 

IVOUT An array of up to 50 elements specifying which variables are to be printed. The 

variables currently available for printing are listed and defined in Table 3-3. 4 The 
default values are 1, 2, 3, 20, 30, 40, 44*0, corresponding to printout of x, y, and 
z- velocity, and static density, pressure, and temperature. 

IWOUT1 A 2-element array, given as IWOUTl(IBOUND), indicating whether or not vari- 
ous parameters are to be printed along the £ boundaries. The subscript 
IBOUND = 1 or 2, corresponding to the £ = 0 and ( = 1 boundaries, respectively. 
Valid values of IWOUTl(IBOUND) are: 

0 for no printout. 

1 for printout along the boundary, with normal derivatives computed using 
three-point one-sided differencing. 

— 1 for printout along the boundary, with normal derivatives computed using 
two-point one-sided differencing. 


4 The definitions of k t and k t in Table 3-3 (IVOUT = 92 and 102) assume a constant turbulent Prandtl number is 
being specified in namelist TURB, 
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The default values are both 0. 


IWOUT2 


IWOUT3 


IDEBUG 


A 2-element array, given as I\VOUT2(IBOUND), indicating whether or not vari- 
ous parameters are to be printed along the rj boundaries. The subscript 
IBOUND = 1 or 2, corresponding to the >7 = 0 and r\ = 1 boundaries, respectively. 
Valid values of IWOUT2(IBOUND) are: 

0 for no printout. 

1 for printout along the boundary, with normal derivatives computed using 
three-point one-sided differencing. 

— 1 for printout along the boundary, with normal derivatives computed using 

two -point one-sided differencing. 

The default values are both 0. 

A 2-element array, given as I\VOUT3(IBOUND), indicating whether or not vari- 
ous parameters are to be printed along the £ boundaries. The subscript 
IBOUND = 1 or 2, corresponding to the £ — 0 and £ = 1 boundaries, respectively. 
Valid values of IWOUT3(IBOUND) are: 

0 for no printout. 

1 for printout along the boundary, with normal derivatives computed using 
three-point one-sided differencing. 

— 1 for printout along the boundary, with normal derivatives computed using 

two-point one-sided differencing. 

The default values are both 0. 

The parameters printed via the IWOUT1, IWOUT2, and/or IWOUT3 input ar- 
rays are the Cartesian coordinates x, y , and z, the static pressure p } the skin friction 
coefficient c h the shear stress r Wi the static temperature T, the heat transfer coeffi- 
cient A, the heat flux q w , and the Stanton number St . See Section 4. 1 for the defi- 
nitions of these parameters. Note that some of these parameters are meaningful 
only if the boundary is a solid wall. In general, using three-point one-sided differ- 
ences to compute normal derivatives will be more accurate. However, for nearly- 
separated flows, it has been found that three-point differencing can give misleading 
negative values of r w and c y. The parameters are printed at the boundary points 
specified by the input parameters IPRT1, IPRT2, and IPRT3, or IPRT1A, 
IPRT2A, and IPRT3A. 

An array of up to 20 elements used to turn on additional printout, normally used 
for debugging purposes. Except where noted, set IDEBUG(I) = 1 for printout 
number I. For options 1 through 7, the input parameters IPRT1, IPRT2, and 
IPRT3, or IPRT1A, IPRT2A, and IPRT3A, determine the grid points at which 
the printout appears. Note that some of these options can generate a lot of output. 
Judicious use of the TPRT" controls is recommended. The debug options cur- 
rently available are as follows: 

Number Printout 

1 Coefficient block submatrices and source term subvectors at time 
level a = IDEBUG(1) if IDEBUG(1)>0, or at time levels 
n> (IDEBUG(l)l if IDEBUG(l) < 0. This printout is done af- 
ter the elimination of any off-diagonal boundary 7 condition sub- 
matrices (subroutine BCELIM) and after any artificial viscosity 
has been added (subroutine AVISC1 or 2), but before any rear- 
rangement of the elements in the boundary condition submatrices 
(subroutine FILTER.) 
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IUNITS 

IPRT 

IPRTA 

IPRT1 

IPRT2 

IPRT3 


Number Printout 

2 Coefficient block submatrices and source term subvectors at time 
level n = IDEBUG(2) if IDEBUG(2)>0, or at time levels 
n > | IDEBUG(2) | if IDEBUG(2) < 0. This printout is done af- 
ter the elimination of any off-diagonal boundary condition sub- 
matrices (subroutine BCELIM), but before any artificial viscosity 
has been added (subroutine AVISC1 or 2) and before any rear- 
rangement of the elements in the boundary' condition submatrices 
(subroutine FILTER.) 

3 Boundary condition coefficient block submatrices and source term 
subvectors at time level 77 = IDEBUG(3) if IDEBUG(3) > 0, or 
at time levels n > | IDEBUG(3) | if IDEBUG(3) < 0. This print- 
out is done before the elimination of any off-diagonal boundary 
condition submatrices (subroutine BCELIM) and before any re- 
arrangement of the elements in the boundary condition submatri- 
ces (subroutine FILTER.) 

4 Boundary condition coefficient block submatrices and source term 
subvectors at time level 77 = IDEBUG(4) if IDEBUG(4) > 0, or 
at time levels 77 > | IDEBUG(4) | if IDEBUG(4) < 0. This print- 
out is done after the elimination of any off-diagonal boundary 
condition submatrices (subroutine BCELIM) and after any rear- 
rangement of the elements in the boundary condition submatrices 
(subroutine FILTER.) 

5 Intermediate solutions Q* and Q** at time level 72 = IDEBUG(5) 
if IDEBUG(5) > 0, or at time levels 72 > |IDEBUG(5)| if 
IDEBUG(5) < 0. 

6 Final solution Q" after the last ADI sweep at time level 
n = IDEBUG(6) if IDEBUG(6) > 0, or at time levels 
77 > | IDEBUG(6) | if IDEBUG(6) < 0. 

7 Cartesian coordinates, metric coefficients, and inverse of the grid 
transformation Jacobian computed in subroutine METS, 

The default values are all 0, 


0 for input and output in English units. 

1 for input and output in SI units. 

The default value is 0. 

Results are printed every IPRT'th time level. However, the initial and final flow 
fields are always printed. The default value is 1. 

An array of up to 101 elements specifying the time levels at which results are to 
be printed. The initial conditions are at time level 1. If the calculation converges, 
or if the pressure or temperature is non-positive, the results are printed regardless 
of the value of IPRTA. If this parameter is specified, it overrides the value of 
IPRT. The default values are all 0. 

Results are printed at every IPRTLth grid point in the £ direction. However, the 
results at the boundaries are always printed. The default value is 1. 

Results are printed at every IPRT2'th grid point in the r\ direction. However, the 
results at the boundaries are always printed. The default value is 1. 

Results are printed at every IPRT3'th grid point in the £ direction. However, the 
results at the boundaries are always printed. The default value is L 


16 3.0 Input Description 


Proteus 3-D User's Guide 



IPRT1A 


An array of up to N1 elements (see Namelist NUM) specifying the { indices at 
which results are to be printed. The indices must be specified in ascending order. 
If this parameter is specified, it overrides the value of IPRT1. The default values 
are all 0. 

IPRT2A An array of up to N2 elements (see Namelist NUM) specifying the r\ indices at 

which results are to be printed. The indices must be specified in ascending order. 
If this parameter is specified, it overrides the value of IPRT2. The default values 
are all 0. 

IPRT3A An array of up to N3 elements (see Namelist NUM) specifying the £ indices at 

which results are to be printed. The indices must be specified in ascending order. 
If this parameter is specified, it overrides the value of IPRT3. The default values 
are all 0. 

NHMAX Maximum number of time levels allowed in the printout of the convergence history 
file (not counting the first two, which are always printed.) The default value is 100. 

Plot File Controls 

In addition to the printed output, files called plot files may be written for use by various post-processing 
routines. The following parameters specify the type of plot files to be written, and at what locations in both 
time and space. These plot files are described in greater detail in Section 4.2. 

IPLOT 0 for no plot file. 

1 to write results into an auxiliary file, in CONTOUR format, for later post- 
processing. If multiple time levels are to be written into the file, they will be 
stacked sequentially. The value of the time t will not be written into the file. 

2 to write results into auxiliary files, in PLOT3D/WHOLE format. Multiple 
time levels will be stacked sequentially, 5 with stored in the Q file 
header, 6 

3 to write results into auxiliary files, in PLOT3D/PLANES format. Multiple 
time levels will be stacked sequentially 5 with Ti tlj i stored in the Q file 
header. 6 

The default value is 0. 

IPLT Results are written into the plot file every IPLT'th time level. However, if 

IPLT > 0, the initial and final flow fields are automatically included in the file. If 
IPLT = 0, only the final flow field is written into the file. The default value is 0. 

IPLT A An array of up to 101 elements specifying the time levels at which results are to 

be written into the plot file. The initial conditions are at time level 1. If the cal- 
culation converges, or if the pressure or temperature is non-positive, the results are 
written into the plot file regardless of the value of IPLT A. If this parameter is 
specified, it overrides the value of IPLT. The default values are all 0. 

Unit Numbers 

The following parameters specify the Fortran unit numbers used for various input and output files. 
NIN, the unit number for reading the standard input file, is hardwired in the program as 5. 

NOUT Unit number for printing standard output. The default value is 6. 


5 The current version of PLOT3D does not work for multiple time levels, although future versions might. 

6 Note that with IDTAU = 5 or 6, t will vary in space, and therefore ^ t^ij. 
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NGRID Unit number for reading computational coordinate system file. The default value 

is 7. 

NPLOTX Unit number for writing XYZ file when using PL0T3D plot file format. The de- 
fault value is 8. 

NPLOT Unit number for writing CONTOUR plot file, or for writing Q file when using 

PLOT3D format. The default value is 9. 

NHIST Unit number for writing convergence history file. The default value is 10. 

3. 1 .5 Namelist GMTRY 


Coordinate System Type 

These parameters specify the type of flow domain being analyzed. Simple Cartesian or cylindrical con- 
figurations can be done automatically. For more complex geometries, the configuration is determined by 
reading a pre-stored coordinate file. Note however, that the number of grid points and their distribution 
can be changed by the parameters in namelist NUM. 

NGEOM Flag used to specify type of computational coordinates. Currently coded are: 

1 Cartesian (x-y-z) computational coordinates. 

2 Cylindrical (r-#-x) computational coordinates. 

10 Get computational coordinates from coordinate system file. The contents of 
this file are described in Section 3.2. 

The default value is 1 . 


Cartesian Computational Coordinates 

The following parameters specify the size of the flow domain for the Cartesian coordinate option 
vNGEOM = 1). The computational (£,»;,£) domain for this option is shown in physical (xj/,z) space in 
Figure 3.1. The points in the figure labeled "min" and "max" are the points (x mini y miny z min ) and 
{Xm CX ,ymcx>Zmc X ), respectively. 

XMIN Minimum x-coordinate for Cartesian coordinate option. The default value is 0.0. 

XMAX Maximum x-coordinate for Cartesian coordinate option. The default value is 1.0. 

YMIN Minimum y-coordinate for Cartesian coordinate option. The default value is 0.0. 

YMAX Maximum y-coordinate for Cartesian coordinate option. The default value is 1.0. 

ZMIN Minimum z-coordinate for Cartesian coordinate option. The default value is 0.0. 

ZMAX Maximum z-coordinate for Cartesian coordinate option. The default value is 1.0. 
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max 



Figure 3.1 - Cartesian computational coordinates. 

Cylindrical Computational Coordinates 

The following parameters specify the size of the flow domain for the cylindrical coordinate option 
(NGEOM = 2). The computational (£,>7,0 domain for this option is shown in physical space in 

Figure 3.2. The points in the figure labeled "min" and "max" are the points (r mm , d mtn , x mm ) and 
(r mcx , 6 max , x max ), respectively. 


RMIN Minimum r-c oordinate for cylindrical coordinate option. The default value is 0.0. 

RMAX Maximum r-coordinate for cylindrical coordinate option. The default value is 1.0. 

THMIN Minimum 0-coordinate in degrees for cylindrical coordinate option. The default 

value is 0.0. 

THMAX Maximum 0-coordinate in degrees for cylindrical coordinate option. The default 
value is 90.0. 

XMIN Minimum ^-coordinate for cylindrical coordinate option. The default value is 0.0. 

XMAX Maximum ^-coordinate for cylindrical coordinate option. The default value is 1.0. 
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Figure 3.2 - Cylindrical computational coordinates. 


3.1.6 Namelist FLOW 


Control^ Flags 

The following parameters are flags that specify the type of equations to be solved, and which variables 
are being supplied as initial conditions. 

IEULER 0 for a full time-averaged Navier-Stokes calculation. 

1 for an Euler calculation (i.e., neglecting all viscous and heat conduction terms.) 

The default value is 0. 

ITHIN A 3-element array, specified as ITHIN(IDIR), indicating whether or not the thin- 

layer option is to be used in direction IDIR. The subscript IDIR = 1, 2, or 3, 
corresponding to the £, tj, and C directions, respectively. Valid values of 
ITHIN(IDIR) are: 

0 to include second derivative viscous terms in direction IDIR. 

1 to use the thin-layer option in direction IDIR. This does not decrease the ex- 
ecution time much, but may be useful if the grid in direction IDIR is not suf- 
ficiently dense to resolve second derivatives in that direction. 


20 3.0 Input Description 


Proteus 3-D User's Guide 



Only one of the three values may be set equal to 1. The default values are all 0. 

IHSTAG 0 to solve the energy equation. The dimensioning parameter XEQP must be 
equal to 5. (See Section 6.2.) 

1 to eliminate the energy equation by assuming constant stagnation enthalpy per 
unit mass. This significantly lowers the overall execution time. 

The default value is 0. 

ILAMV 0 for constant laminar viscosity and thermal conductivity coefficients equal to 

MUR and KTR. 

1 for variable laminar viscosity and thermal conductivity coefficients, computed 
as a function of local temperature using Sutherland's formula for air (White, 
1974). 

The default value is 0. 

ICVARS Parameter specifying which variables are being supplied as initial conditions for the 

time marching procedure by subroutine INIT. Remember that the initial condi- 
tions must be nondimensionalized by the reference conditions listed in Table 3-2. 
(See Section 5.0 for details on defining initial conditions.) When the energy 
equation is being solved (IHSTAG = 0), the allowed values are: 

ICVARS Variables Supplied By INIT 

1 p, pw, pv, pw, Et 

2 p, u, v, w, T 

3 p, u, v, w, T 

4 p, u, v, w, p 

5 c pt u , v, w, T 

6 p, M, a„, a w , T 

When constant stagnation enthalpy is assumed (IHSTAG = 1), the allowed values 
are: 


ICVARS Variables Supplied By INIT 

1 P, pw, pv, pw 

2 p, u, v, w 

3 p, u, v, w 

5 c p , u y v, w 

6 p t M, a„, a w 

In the above tables, c p , a v , and a* represent static pressure coefficient, flow angle 
in degrees in the x-y plane, and flow angle in degrees in the x-z plane, respectively. 
The default value is 2. 


Reference Conditions 

The following parameters specify the six basic reference conditions for length, velocity, temperature, 
density, viscosity, and thermal conductivity. These reference conditions are used, along with some addi- 
tional reference conditions derived from them, as the nondimensionalizing factors for nondimensional input 
and output parameters. The dimensional reference conditions may be read in using either English or SI 
units, depending on the value of IUNITS. 

LR Reference length Lr in feet (meters). This is a real variable. The default value is 

1 . 0 . 
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UR Reference velocity u r in ft/sec (m/sec). Either UR or MACHR may be specified, 

but not both. The unspecified one will be computed from the remaining reference 
conditions. The default value is a r = (y r R 7 r ) i/2 , the speed of sound at the reference 
temperature. 

MACHR Reference Mach number, M r — u r l(y r RT r ) 112 . Either MACHR or UR may be 
specified, but not both. The unspecified one will be computed from the remaining 
reference conditions. This is a real variable. The default value is 0.0. 


TR 

RHOR 


MUR 


RER 

KTR 


PRLR 


£ 7 uid Properties 


Reference temperature T r in °R (K). The default value is 519.0 °R (288.333 K). 

Reference density p r in lb m /ft 3 (kg/m 3 ). The default value is 0.07645 lb m /ft 3 
(1.22461 kg/m 3 ). 

Reference viscosity p r in lb m /ft-sec (kg/m-sec). Either MUR or RER may be 
specified, but not both. The unspecified one will be computed from the remaining 
reference conditions. This is a real variable. The default value is the viscosity for 
air at the reference temperature TR, 

Reference Reynolds number, Re r = p r u r Lrlp r . Either RER or MUR may be spec- 
ified, but not both. The unspecified one will be computed from the remaining 
reference conditions. The default value is 0.0. 

Reference thermal conductivity k r in lb m -ft/sec 3 -°R (kg-m/sec 3 -K). Either KTR 
or PRLR may be specified, but not both. The unspecified one will be computed 
from the remaining reference conditions. This is a real variable. The default value 
is the thermal conductivity for air at the reference temperature TR. 

Reference laminar Prandtl number Pr lr = Cp^lK^ Here c Pr is the specific heat at 
constant pressure, defined as either c p (T r ) or y r RI(y r — 1), depending on whether 
or not a value is being specified for the input parameter GAMR. Either PRLR 
or KTR may be specified, but not both. The unspecified one will be computed 
from the remaining reference conditions. The default value is 0.0. 


The following parameters provide information about the fluid being used. 

RG Gas constant R in ft 2 /sec 2 -°R (m 2 /sec 2 -K). The default value is 1716 ft 2 /sec 2 -°R 

(286.96 m 2 /sec 2 -K). 

GAMR Reference ratio of specific heats, y r — c Pr jc Vr . This parameter acts as a flag for a 

constant specific heat option. If a non-zero value for GAMR is specified by the 
user, the specific heats c p and c v are computed from GAMR and RG, and treated 
as constants. Otherwise they are computed locally as a function of temperature. 
The default value is 0.0. 


HSTAGR Stagnation enthalpy h T in ft 2 /sec 2 (m 2 /sec 2 ). This parameter is only used with the 
constant stagnation enthalpy option (IHSTAG =1). The default value is com- 
puted from the reference conditions. 

3.1.7 Namelist BC 


The parameters in this namelist specify the boundary conditions to be used for the mean flow equations 
and, if necessary, for the k-s turbulence model equations. For the mean flow T , NEQ conditions must be 
specified at each computational boundary, where NEQ is the number of coupled equations being solved. 
NEQ will be equal to 4 or 5 depending on the value of IHSTAG. (See Table 3-4.) 
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Note that the boundary conditions may be thought of as simply NEQ additional equations to be solved 
on the boundary". They do not necessarily have to be associated one-to-one with the governing differential 
equations or the dependent variables. They must, however, be functions of the dependent variables and 
sufficiently complete to set constraints on each of the dependent variables through their functional form. 
They must also, of course, be independent of one another and physically appropriate for the problem being 
solved. 

Three different methods are available for setting boundary conditions for steady flow computations. 
The first, and easiest, way is to specify the type of boundary (i.e., solid wall, symmetry, etc.) using the KBC 
input parameters. These parameters act as Tneta" flags, triggering the automatic setting of the necessary 
NEQ individual boundary' conditions at the specified boundary. 

Second, if more flexibility is needed, the NEQ individual boundary conditions may be set for each 
boundary using the JBC and GBC input parameters. The boundary condition type (specified value, speci- 
fied gradient, etc.) is given by JBC, and the boundary condition value by GBC. With these parameters, the 
same conditions are applied over the entire surface. 

And third, if even greater flexibility is needed, the NEQ individual boundary conditions may be set for 
each boundary using the 1BC and FBC input parameters. These are analogous to the JBC and GBC pa- 
rameters (i.e., the boundary condition type is given by IBC, and the value by FBC), but they allow a 
point-by-point distribution of type and value to be specified instead of using the same type and value over 
the entire surface. 7 

For a given boundary, boundary conditions specified via the KBC parameters override those specified 
using the JBC and GBC parameters, which in turn override those specified using the IBC and FBC pa- 
rameters. However, the different methods may be used in combination as long as they don't conflict. For 
example, the KBC parameters may be used for four boundaries, the JBC and GBC parameters for the fifth 
boundary, and the IBC and FBC parameters for the sixth boundary. And, on a single boundary, the JBC 
and GBC parameters may be used for some of the NEQ boundary conditions, and the IBC and FBC pa- 
rameters for the rest. 

Unsteady boundary conditions may be used w r hen individual boundary conditions are specified for the 
entire surface, but not when boundary conditions are specified point-by-poim.. 

With one exception, the NEQ boundary conditions at each boundary may be specified in any order. 
The exception is any condition on one of the dependent conservation variables Q. These must be specified 
in the order given in Table 3-4. 

If a problem requires a boundary condition of the form A F = 0 , F=f dFjd4> =/, or VF- n =/, where 
F is not one of the functions already built into Proteus , the subroutines BCF and BCFLIN may be used. 
This requires that the user supply subroutine BCFLIN. Subroutines BCF and BCFLIN are described in 
detail in Volume 3. 

If the Chien k-z turbulence model is being used, boundary 7 conditions for k and z must be specified in 
addition to those specified for the mean flow equations. (Unless a spatially periodic boundary condition is 
being used for the mean flow. In this case, separate boundary conditions for k and z are not needed.) 
Surface boundary conditions may be specified using the JBCT/GBCT parameters, or point-by-point 
boundary conditions may be specified using the IBCT/FBCT parameters. These parameters are analogous 
to the JBC/GBC and IBC/FBC parameters used for the mean flow boundary conditions. There are no 
input parameters for k-z boundary conditions that are analogous to the KBC parameters. 


7 However, note that a specified point-by-point distribution of a function value is most easily set using the 'no change 
from initial conditions' option with the JBC parameters. 
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Mean Flow Boundary Types with KBC 


The following parameters set mean flow boundary conditions by specifying the type of boundary (i.e., 
solid wall, symmetry, etc.)- These parameters act as "meta" flags, triggering the automatic setting of the 
necessary JBC and GBC values. 

KBC1 An array, given as KBCl(IBOUXD), specifying the types of boundaries in the £ 

direction. The subscript IBOUND = I or 2, corresponding to the £ = 0 and 
£ = 1 boundaries, respectively. The default values are both 0. 

KBC2 An array, given as KBC2(IBOUND), specifying the types of boundaries in the n 

direction. The subscript IBOUND = 1 or 2, corresponding to the tj = 0 and 
jj = 1 boundaries, respectively. The default values are both 0. 

KBC3 An array, given as KBC3(IBOUND), specifying the types of boundaries in the £ 

direction. The subscript IBOUND =1 or 2, corresponding to the £ = 0 and 
£ = 1 boundaries, respectively. The default values are both 0. 

The boundary' types that may be specified are described briefly in the following table, and in greater 
detail in Table 3-5. For boundary types involving gradient boundary conditions, 2-point differencing is used 
if the input KBC value is positive, and 3-point differencing is used if it is negative. For boundary' types 
involving "no change from initial conditions"-type boundary conditions (e.g., AT = 0), the proper boundary- 
values must be set in the initial conditions. When the KBC parameters are used, the corresponding GBC 
and FBC parameters should be defaulted. 


KBC Value Boundary Type 


± 1 

No-slip adiabatic wall- 

+ 2 

No-slip wall, specified temperature. 

± 3 

Inviscid wall. 

10 

Subsonic inflow, linear extrapolation. 

+ 11 

Subsonic inflow, zero gradient. 

20 

Subsonic outflow, linear extrapolation. 

+ 21 

Subsonic outflow, zero gradient. 

30 

Supersonic inflow. 

40 

Supersonic outflow, linear extrapolation. 

+ 41 

Supersonic outflow, zero gradient. 

+ 50 

Symmetry. 

60 

Spatially periodic. 


Boundary conditions specified using the KBC parameter for a given boundary override any boundary- 
condition types specified for that boundary using the JBC or IBC parameters. Note, however, that since 
the default values for the KBC parameters are all 0, the default procedure for specifying boundary conditions 
is by using the JBC and GBC parameters. 

Surface Mean Flow Bomdari Condition T wes and_ Values^ with JBC and GBC 

The following parameters set the NEQ individual mean flow boundary condition types and values for 
each boundary using the JBC and GBC parameters. With these parameters, the same conditions are applied 
over the entire surface. Remember that the boundary condition values must be nondimensionalized by the 
reference conditions listed in Table 3-2. If some of the boundary conditions are being specified using the 
IBC and FBC parameters, the appropriate JBC parameters must be set equal to — 1, as described below. 

JBC1 A two-dimensional array, given as JBC 1(IEQ, IBOUND), specifying the type of 

boundary' conditions to be used on the £ = 0 and £ = 1 boundaries. Here 
IEQ = 1 to NEQ corresponding to each equation, and IBOUND = 1 or 2 corre- 
sponding to the £ = 0 and £ = 1 boundaries, respectively. Setting JBC1 = — 1 
signals the code to use boundary conditions specified point-by-point, as given by 
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the input arrays IBC1 and FBC1. See Table 3-6 for a list of allowed boundary 
condition types. The default values are all 0. 


JBC2 A two-dimensional array, given as JBC2(IEQ,IBOUND), specifying the type of 

boundary conditions to be used on the ly — 0 and yj = 1 boundaries. Here 
IEQ = 1 to NEQ corresponding to each equation, and IBOUND = 1 or 2 corre- 
sponding to the v\ = 0 and rj = 1 boundaries, respectively. Setting JBC2 = - 1 
signals the code to use boundary conditions specified point-by-point, as given by 
the input arrays IBC2 and FBC2. See Table 3-6 for a list of allowed boundary 
condition types. The default values are all 0. 

JBC3 A two-dimensional array, given as JBC3(IEQ, IBOUND), specifying the type of 

boundary conditions to be used on the £ = 0 and £ = 1 boundaries. Here 
IEQ = 1 to NEQ corresponding to each equation, and IBOUND = 1 or 2 corre- 
sponding to the £ — 0 and £ = 1 boundaries, respectively. Setting JBC3 = - 1 
signals the code to use boundary conditions specified point-by-point, as given by 
the input arrays IBC3 and FBC3. See Table 3-6 for a list of allowed boundary 
condition types. The default values are all 0. 

GBC1 A two-dimensional array, given as GBC 1 (IEQ, IBOUND), specifying the values for 

the steady boundary conditions to be used on the £ = 0 and £ = 1 boundaries. 
Here IEQ = 1 to NEQ corresponding to each equation, and IBOUND = 1 or 2 
corresponding to the £ = 0 and £ = 1 boundaries, respectively. The default values 
are all 0 . 0 . 

GBC2 A two-dimensional array, given as GBC2(IEQ, IBOUND), specifying the values for 

the steady boundary conditions to be used on the rj = 0 and > 7=1 boundaries. 
Here IEQ = 1 to NEQ corresponding to each equation, and IBOUND = 1 or 2 
corresponding to the >7 = 0 and > 7=1 boundaries, respectively. The default values 
are all 0 . 0 . 

GBC3 A two-dimensional array, given as GBC 3 (IEQ, IBOUND), specifying the values for 

the steady boundary conditions to be used on the £ = 0 and £ = l boundaries. 
Here IEQ = 1 to NEQ corresponding to each equation, and IBOUND = I or 2 
corresponding to the £ = 0 and £ = 1 boundaries, respectively. The default values 
are all 0 . 0 . 

Note that boundary condition types 2, 12, 22, etc., are specified values of the derivative with respect to 
the computational coordinate, not with respect to the physical distance in the direction of the computational 
coordinate. See Section 6.3 of Volume 1 for details. 

Note also that normal derivative boundary condition values are positive in the direction of increasing 
£ rj, or £. Thus, a positive value for GBC at £ = 0, >7 = 0, or £ = 0 implies a flux into the computational 
domain, and a positive GBC at £ = 1, yj = 1, or £ = 1 implies a flux out of the computational domain. See 
Section 6.4 of Volume 1 for details. Similarly, the normal velocity V n is positive in the direction of m- 
creasing £, >7, or £. Thus, a positive V n at £ = 0, rj = 0, or £ = 0 implies flow into the computational domain, 
and a positive V n at £ = 1, 17 = 1, or £ = 1 implies flow out of the computational domain. See the de- 
scription of subroutine BCVDIR in Section 4.2 of Volume 3 for details. 


Boundary conditions specified using the JBC and GBC parameters for given values of IEQ and 
IBOUND override any boundary conditions specified for those values of IEQ and IBOUND using the IBC 
and FBC parameters. Note that since the default values for the JBC parameters are all 0, the default 
boundary conditions are "no change from initial conditions" for the conservation variables. 

Point-by- Point Mean Flow Boundary Condition Types and Values with IBC and FBC 

The following parameters set the NEQ individual mean flow boundary condition types and values for 
each boundary using the IBC and FBC parameters. With these parameters, point-by-point distributions 
are specified on the surface for the boundary condition types and values. Remember that the boundary 
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condition values must be nondimensionalized by the reference conditions listed in Table 3-2. Note that 
these parameters are activated by setting the appropriate JBC parameters equal to — 1, as described below. 

IBC1 A four-dimensional array, given as IBCl(I2,I3,IEQ,IBOUND), specifying the type 

of boundary conditions to be used on the £ = 0 and £ = 1 boundaries. Here 
12 = 1 to N2 and 13 = 1 to N3 corresponding to each grid point on the boundary, 
IEQ = I to NEQ corresponding to each equation, and IBOUND = 1 or 2 corre- 
sponding to the £ = 0 and £ = 1 boundaries, respectively. JBC 1(IEQ, IBOUND) 
must be set equal to — 1. See Table 3-6 for a list of allowed boundary condition 
types. The default values are all 0. 

IBC2 A four-dimensional array, given as IBC2(1 1,13, IEQ, IBOUND), specifying the type 

of boundary conditions to be used on the y\ = 0 and r\ = 1 boundaries. Here 
II = 1 to N1 and 13 = 1 to N3 corresponding to each grid point on the boundary, 
IEQ = 1 to NEQ corresponding to each equation, and IBOUND = 1 or 2 corre- 
sponding to the rj — 0 and r\ = 1 boundaries, respectively. JBC2(IEQ, IBOUND) 
must be set equal to - 1. See Table 3-6 for a list of allowed boundary condition 
types. The default values are all 0. 

IBC3 A four-dimensional array, given as IBC3(I1, 12, IEQ, IBOUND), specifying the type 

of boundary conditions to be used on the £ = 0 and £ = 1 boundaries. Here 
II * 1 to N1 and 12 = 1 to N2 corresponding to each grid point on the boundary, 
IEQ = 1 to NEQ corresponding to each equation, and IBOUND = 1 or 2 corre- 
sponding to the £ = 0 and £ = 1 boundaries, respectively. JBC3(IEQ, IBOUND) 
must be set equal to - 1. See Table 3-6 for a list of allowed boundary condition 
types. The default values are all 0. 

FBC1 A four-dimensional array, given as FBC1(I2, 13, IEQ, IBOUND), specifying the 

values for the steady boundary conditions to be used on the £ = 0 and £ = 1 
boundaries. Here 12 = 1 to N2 and 13 = 1 to N3 corresponding to each grid point 
on the boundary, IEQ = 1 to NEQ corresponding to each equation, and 
IBOUND = 1 or 2 corresponding to the £ = 0 and £ = 1 boundaries, respectively. 
The default values are all 0.0. 

FBC2 A four-dimensional array, given as FBC2(I1, 13, IEQ, IBOUND), specifying the 

values for the steady boundary conditions to be used n the r\ — 0 and tj — 1 
boundaries. Here II = 1 to N1 and 13 = 1 to N3 corresponding to each grid point 
on the boundary, IEQ=1 to NEQ corresponding to each equation, and 
IBOUND = 1 or 2 corresponding to the q = 0 and tj = 1 boundaries, respectively. 
The default values are all 0.0. 

FBC3 A four-dimensional array, given as FBC3(I1, 12, IEQ, IBOUND), specifying the 

values for the steady boundary conditions to be used on the £ = 0 and £ = 1 
boundaries. Here II = 1 to N1 and 12 = 1 to N2 corresponding to each grid point 
on the boundary, IEQ = 1 to NEQ corresponding to each equation, and 
IBOUND = 1 or 2 corresponding to the £ = 0 and £ = 1 boundaries, respectively. 
The default values are all 0.0. 

Note that boundary condition types 2, 12, 22, etc., are specified values of the derivative with respect to 
the computational coordinate, not with respect to the physical distance in the direction of the computational 
coordinate. See Section 6.3 of Volume 1 for details. 

Note also that normal derivative boundary condition values are positive in the direction of increasing 
£, rj, or £. Thus, a positive value for FBC at £ = 0, rj = 0, or £ = 0 implies a flux into the computational 
domain, and a positive FBC at £ = 1, y\ — 1, or £ = 1 implies a flux out of the computational domain. See 
Section 6.4 of Volume 1 for details. Similarly, the normal velocity V n is positive in the direction of in- 
creasing £, rj, or £. Thus, a positive V„ at £ = 0, rj = 0, or £ = 0 implies flow into the computational domain, 
and a positive V„ at £ = 1, rj = 1, or £ = 1 implies flow out of the computational domain. See the de- 
scription of subroutine BCVDIR in Section 4.2 of Volume 3 for details. 
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Unsteady Mean Flow » Boundary Conditions 


The following parameters are used to specify unsteady mean flow boundary conditions. The boundary 
condition type (specified value, specified gradient, etc.) is given by JBC, as described above, but the value 
is given by GTBC. The type of unsteadiness (general or periodic) is given by JTBC. 

JTBC1 A two-dimensional array, given as JTBCl(IEQ,IBOUND), specifying the type of 

time dependency for the boundary conditions on the c = 0 and q = 1 boundaries. 
Here IEQ = 1 to NEQ corresponding to each equation, and IBOUND = 1 or 2 
corresponding to the f = 0 and £ = 1 boundaries, respectively. Valid values of 
JTBC 1(IEQ, IBOUND) are: 

0 for a steady boundary condition, whose value is given by GBC1. 

1 for a general unsteady boundary condition, whose value is determined by linear 
interpolation in the input table of GTBC1 vs. NTBCA. 

2 for a time-periodic boundary condition of the form g\ + g 2 sin(g3/? 4- g 4), where 
n is the time level and g { through g 4 are given by the first four values of GTBC 1 . 

The default values are all 0. 

JTBC2 A two-dimensional array, given as JTBC2(IEQ, IBOUND), specifying the type of 

time dependency for the boundary conditions on the tj - 0 and rj = 1 boundaries. 
Here IEQ = 1 to NEQ corresponding to each equation, and IBOUND = 1 or 2 
corresponding to the rj = 0 and rj = 1 boundaries, respectively. Valid values of 
JTBC2(IEQ, IBOUND) are: 

0 for a steady boundary condition, whose value is given by GBC2. 

1 for a general unsteady boundary condition, whose value is determined by linear 
interpolation in the input table of GTBC2 vs. NTBCA. 

2 for a time-periodic boundary condition of the form g\ 4- gi sin (g^n 4* g 4 ), where 
n is the time level and g x through g 4 are given by the first four values of GTBC2. 

The default values are all 0. 

JTBC3 A two-dimensional array, given as JTBC3(IEQ, IBOUND), specifying the type of 

time dependency for the boundary conditions on th, £ = 0 and C = 1 boundaries. 
Here IEQ = 1 to NEQ corresponding to each equation, and IBOUND = 1 or 2 
corresponding to the £ = 0 and £ = 1 boundaries, respectively. Valid values of 
JTBC3(IEQ, IBOUND) are: 

0 for a steady boundary condition, whose value is given by GBC3. 

1 for a general unsteady boundary condition, whose value is determined by linear 
interpolation in the input table of GTBC3 vs. NTBCA. 

2 for a time-periodic boundary condition of the form g\ 4- gi sin (g^n 4- g 4 ), where 
n is the time level and g x through g 4 are given by the first four values of GTBC3. 

The default values are all 0. 

NTBC Number of values in the tables of GTBC1, GTBC2, and/or GTBC3 vs. NTBCA 

for the general unsteady boundary condition option. The maximum value allowed 
is the value of the dimensioning parameter NTP. (See Section 6.2.) The default 
value is 0. 

NTBCA An array of NTBC time levels at which GTBC1, GTBC2, and/or GTBC3 are 

specified for the general unsteady boundary condition option. The default values 
are all 0. 

GTBC1 A three-dimensional array, given as GTBC 1(ITBC, IEQ, IBOUND), used in the 

unsteady and time -periodic boundary condition options for the £ = 0 and £ = 1 
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boundaries. Here IEQ = 1 to XEQ corresponding to each equation, and 
IBOUND = 1 or 2 corresponding to the £ = 0 and £ = 1 boundaries, respectively. 
For general unsteady boundary conditions the subscript ITBC = 1 to NTBC, 
corresponding to the time levels in the array NTBCA, and GTBC1 specifies the 
boundary condition value directly. For time-periodic boundary 7 conditions the 
subscript ITBC — 1 to 4, and GTBC1 specifies the four coefficients in the equation 
used to determine the boundary condition value. The default values are all 0.0. 

GTBC2 A three-dimensional array, given as GTBC2(ITBC,IEQ J IBOUND) ) used in the 

unsteady and time-periodic boundary condition options for the rj = 0 and y[ = 1 
boundaries. Here IEQ = 1 to NEQ corresponding to each equation, and 

IBOUND = 1 or 2 corresponding to the >7 = 0 and r\ = 1 boundaries, respectively. 
For general unsteady boundary conditions the subscript ITBC = 1 to NTBC, 
corresponding to the time levels in the array NTBCA, and GTBC2 specifies the 
boundary condition value directly. For time-periodic boundary conditions the 
subscript ITBC = 1 to 4, and GTBC2 specifies the four coefficients in the equation 
used to determine the boundary 7 condition value. The default values are all 0.0. 

GTBC3 A three-dimensional array, given as GTBC3(ITBC,IEQ,IBOUND), used in the 

unsteady and time -periodic boundary condition options for the £ = 0 and £ = 1 
boundaries. Here IEQ = 1 to NEQ corresponding to each equation, and 

IBOUND = 1 or 2 corresponding to the £ = 0 and £ = 1 boundaries, respectively. 
For general unsteady boundary conditions the subscript ITBC = 1 to NTBC, 
corresponding to the time levels in the array NTBCA, and GTBC3 specifies the 
boundary condition value directly. For time-periodic boundary conditions the 
subscript ITBC = 1 to 4, and GTBC3 specifies the four coefficients in the equation 
used to determine the boundary condition value. The default values are all 0.0. 

k-z Surface Boundary Condition Types and Values with JBCT and GBCT 

The following parameters set the individual k and e boundary condition types and values for each 
boundary using the JBCT and GBCT parameters. With these parameters, the same conditions are applied 
over the entire boundary. Remember that the boundary condition values must be nondimensionalized by 
the reference conditions listed in Table 3-2. None of the following parameters are needed if spatially peri- 
odic boundary conditions are being used. If either of the boundary conditions for a computational 
boundary is being specified using the IBCT and FBCT parameters, the appropria 1 JBCT parameters must 
be set equal to — 1, as described below, 

JBCT1 A two-dimensional array, given as JBCT1(IEQ, IBOUND), specifying the type of 

boundary conditions to be used on the £ = 0 and £ = 1 boundaries. Here 
IEQ = 1 or 2 corresponding to k and e, respectively, and IBOUND = 1 or 2 cor- 
responding to the £ — 0 and £ = 1 boundaries, respectively. Setting JBCT1 = — 1 
signals the code to use boundary 7 conditions specified point -by- point, as given by 
the input arrays IBCT1 and FBCT1. See Table 3-7 for a list of allowed boundary 
condition types. The default values are 0 for IEQ = 1, and 10 for IEQ = 2. 

JBCT2 A two-dimensional array, given as JBCT2(IEQ, IBOUND), specifying the type of 

boundary conditions to be used on the rj = 0 and *7 = 1 boundaries. Here 
IEQ = 1 or 2 corresponding to k and s, respectively, and IBOUND = 1 or 2 cor- 
responding to the rj = 0 and y\ = 1 boundaries, respectively. Setting JBCT2 = — 1 
signals the code to use boundary conditions specified point-by-point, as given by 
the input arrays IBCT2 and FBCT2. See Table 3-7 for a list of allowed boundary 
condition types. The default values are 0 for IEQ = 1, and 10 for IEQ = 2. 

JBCT3 A two-dimensional array, given as JBCT3(IEQ, IBOUND), specifying the type of 

boundary 7 conditions to be used on the £ = 0 and £ = 1 boundaries. Here 
IEQ = 1 or 2 corresponding to k and e, respectively, and IBOUND = 1 or 2 cor- 
responding to the £ — 0 and £ = 1 boundaries, respectively. Setting JBCT3 = — 1 
signals the code to use boundary conditions specified point-by-point, as given by 
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the input arrays IBCT3 and FBCT3. See Table 3-7 for a list of allowed boundary 
condition types. The default values are 0 for IEQ = 1, and 10 for IEQ = 2. 

GBCT1 A two-dimensional array, given as GBCT 1 (IEQ, IBOUND), specifying the values 

for the boundary conditions to be used on the £ = 0 and £ = 1 boundaries. Here 
IEQ = 1 or 2 corresponding to k and z, respectively, rind IBOUND = 1 or 2 cor- 
responding to the £ = 0 and £ = 1 boundaries, respectively. The default values are 
all 0 . 0 . 

GBCT2 A two-dimensional array, given as GBCT2(IEQ, IBOUND), specifying the values 

for the boundary conditions to be used on the y = 0 and 77 = 1 boundaries. Here 
IEQ = 1 or 2 corresponding to k and e, respectively, and IBOUND = 1 or 2 cor- 
responding to the rj = 0 and >7 = 1 boundaries, respectively. The default values are 
all 0 . 0 . 

GBCT3 A two-dimensional array, given as GBCT3(IEQ, IBOUND), specifying the values 

for the boundary conditions to be used on the £ = 0 and £ = 1 boundaries. Here 
IEQ = 1 or 2 corresponding to k and z, respectively, and IBOUND = 1 or 2 cor- 
responding to the £ = 0 and £ = 1 boundaries, respectively. The default values are 
all 0 . 0 . 

Note that boundary condition types 2 and 12 are specified values of the derivative with respect to the 
computational coordinate, not with respect to the physical distance in the direction of the computational 
coordinate. See Section 6.3 of Volume 1 for details. 

Boundary conditions specified using the JBCT and GBCT parameters for given values of IEQ and 
IBOUND will override any boundary' conditions specified for those values of IEQ and IBOUND using the 
IBCT and FBCT parameters. Note that the default values for the JBCT parameters are "no change from 
initial conditions" for the k and z. 

k-z Point-by- Point Boundary Condition Tyner and Valuer with IBCT and FBCT 

The following parameters set the individual k and t boundary condition types and values for each 
boundary using the IBCT and FBCT parameters. With these parameters, point-by-point distributions are 
specified on the boundary for the boundary condition types and values. Remember that the boundary 
condition values must be nondimensionalized by the reference conditions list, i in Table 3-2. None of the 
following parameters are needed if spatially periodic boundary conditions are being used. Note that these 
parameters are activated by setting the appropriate JBCT parameters equal to - I, as described below. 

IBCT1 A four-dimensional array, given as IBCT1(I2,I3,IEQ,IB0UND), specifying the 

type of boundary conditions to be used on the £ = 0 and £ = 1 boundaries. Here 
12 = 1 to N2 and 13 = 1 to N3 corresponding to each grid point on the boundary, 
IEQ = 1 or 2 corresponding to k and z, respectively, and IBOUND = 1 or 2 cor- 
responding to the £ = 0 and £ = 1 boundaries, respectively. 
JBCT1(IEQ, IBOUND) must be set equal to - 1. See Table 3-7 for a list of al- 
lowed boundary condition types. The default values are 0 for IEQ = 1, and 10 for 
IEQ = 2. 

IBCT2 A four-dimensional array, given as IBCT2(I1, 13, IEQ, IBOUND), specifying the 

type of boundary conditions to be used on the »j = 0 and r; = 1 boundaries. Here 
11 = 1 to N 1 and 13 = 1 to N3 corresponding to each grid point on the boundary, 
IEQ = 1 or 2 corresponding to k and z, respectively, and IBOUND = 1 or 2 cor- 
responding to the y] = 0 and rj = 1 boundaries, respectively. 
JBCT2(IEQ, IBOUND) must be set equal to - 1. See Table 3-7 for a list of al- 
lowed boundary condition types. The default values are 0 for IEQ = 1 , and 1 0 for 
IEQ = 2. 

IBCT3 A four-dimensional array, given as IBCT3(1 1,12, IEQ, IBOUND), specifying the 

type of boundary conditions to be used on the £ = 0 and £ = 1 boundaries. Here 
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II = 1 to XI and 12 = I to X2 corresponding to each grid point on the boundary, 
IEQ = 1 or 2 corresponding to k and s, respectively, and IBOUND = 1 or 2 cor- 
responding to the £ = 0 and £ = I boundaries, respectively. 
JBCT3(IEQ, IBOUND) must be set equal to - 1. See Table 3-7 for a list of al- 
lowed boundary condition types. The default values are 0 for IEQ = 1, and 10 for 
IEQ = 2. 

FBCT1 A four-dimensional array, given as FBCT1(I2, 13, IEQ, IBOUND), specifying the 

values for the boundary conditions to be used on the £ = 0 and £ = 1 boundaries. 
Here 12= 1 to N2 and 13 = 1 to N3 corresponding to each grid point on the 
boundary, IEQ = 1 or 2 corresponding to k and e, respectively, and 

IBOUND = 1 or 2 corresponding to the £ = 0 and £ = 1 boundaries, respectively. 
The default values are all 0.0. 

FBCT2 A four-dimensional array, given as FBCT2(1 1,13, IEQ, IBOUND), specifying the 

values for the boundary conditions to be used on the rj = 0 and >j = 1 boundaries. 
Here 11=1 to N 1 and 13=1 to N3 corresponding to each grid point on the 
boundary, IEQ = 1 or 2 corresponding to k and e, respectively, and 

IBOUND = 1 or 2 corresponding to the ij = 0 and = 1 boundaries, respectively. 
The default values are all 0.0. 

FBCT3 A four-dimensional array, given as FBCT3(I1, 12, IEQ, IBOUND), specifying the 

values for the boundary conditions to be used on the £ = 0 and £ = 1 boundaries. 
Here II = 1 to N1 and 12= 1 to N2 corresponding to each grid point on the 
boundary, IEQ = 1 or 2 corresponding to k and e, respectively, and 

IBOUND = 1 or 2 corresponding to the £ = 0 and £ = 1 boundaries, respectively. 
The default values are all 0.0. 

Note that boundary condition types 2 and 12 are specified values of the derivative with respect to the 
computational coordinate, not with respect to the physical distance in the direction of the computational 
coordinate. See Section 6.3 of Volume 1 for details. 

3.1.8 Namelist N'UM 

Mesh Parameters 

The following parameters specify the number of mesh points and the degree of packing. 

Ml Number of grid points A) in the £ direction. For non-periodic boundary condi- 

tions in the £ direction, the maximum value allowed is the value of the dimen- 
sioning parameter N 1 P. For spatially periodic boundary conditions, the maximum 
is N 1 P — 1 . (See Section 6.2.) The default value is 5. 

N2 Number of grid points N 2 in the >j direction. For non-periodic boundary condi- 

tions in the >j direction, the maximum value allowed is the value of the dimen- 
sioning parameter N2P. For spatially periodic boundary conditions, the maximum 
isN2P— 1. (See Section 6.2.) The default value is 5. 

N3 Number of grid points V 3 in the £ direction. For non-periodic boundary condi- 

tions in the £ direction, the maximum value allowed is the value of the dimen- 
sioning parameter N3P. For spatially periodic boundary conditions, the maximum 
is N3P - 1. (See Section 6.2.) The default value is 5. 

IPACK A 3-element array, specified as IPACK(IDIR), indicating whether or not grid 

points are to be packed in direction IDIR. The subscript IDIR = 1, 2, or 3, 
corresponding to the £, >j, and £ directions, respectively. Valid values of 
IPACK(IDIR) are: 

0 for no packing in direction IDIR. 
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1 to pack points in direction IDIR using a transformation due to Roberts (1971). 
The location and amount of packing are specified by the array SQ. 

The default values are all 0. 

SQ A two-dimensional array controlling the packing of grid points near computational 

boundaries, specified as SQ(IDIR,IPC). The subscript IDIR = 1, 2, or 3 corre- 
sponding to packing in the <*, y\, and £ directions, respectively. The subscript 
IPC = I or 2, where SQ(IDIR,1) specifies the packing location, and SQ(IDIR,2) 
specifies the amount of packing. 

Meaningful values for SQ(IDIR,1) are 0.0, 0.5, and 1.0, where 0.0 corresponds to 
packing near the lower boundary only (i.e., at rj, or £ = 0, depending on IDIR), 
1.0 corresponds to packing near the upper boundary only, and 0.5 corresponds to 
equal packing at both boundaries. 

Meaningful values for SQ(IDIR,2) are values above 1.0, but generally 1. 1 or below. 
The closer SQ(IDIR,2) is to 1.0, the tighter the packing will be. 

The default values are SQ(IDIR,1) = 0.0 and SQ(IDIR,2) = 1.01 for IDIR =1,2, 
and 3. 

Artificial Viscosity Parameters 

The following parameters specify the type and amount of artificial viscosity to be used. 

IAV4E 0 for no fourth-order explicit artificial viscosity. 

1 to include fourth-order explicit artificial viscosity using the constant coefficient 
model of Steger (1978). 

2 to include fourth-order explicit artificial viscosity using the nonlinear coefficient 
model of Jameson, Schmidt, and Turkel (1981). 

The default value is 1 . 

IAV2E 0 for no second-order explicit artificial viscosity. 

1 to include second-order explicit artificial viscosity using the constant coefficient 
model. 

2 to include second-order explicit artificial viscosity using the nonlinear coefficient 
model of Jameson, Schmidt, and Turkel (1981). 

The default value is 0. 

IAV2I 0 for no second-order implicit artificial viscosity. 

1 to include second-order implicit artificial viscosity using the constant coefficient 
model of Steger (1978). 

The default value is 1 . 

CAVS4E For the constant coefficient model, CAVS4E(IEQ) specifies the fourth-order arti- 

ficial viscosity coefficient directly. For the nonlinear coefficient model it speci- 
fies the constant k 4 . The subscript IEQ varies from 1 to NEQ corresponding to 
each coupled equation. (See Table 3-4 for the order of the equations being solved.) 
Good values for a given application are usually determined by experience, but re- 
commended starting values are 1.0 for the constant coefficient model, 0.005 for the 
nonlinear model when spatially varying second-order time differencing is used, and 
0.0002 for the nonlinear model when a spatially constant first-order time differ- 
encing is used. The default values are all 1.0. 
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CAVS2E For the constant coefficient model, CAVS2E(IEQ) specifies the second-order arti- 

ficial viscosity coefficient eg> directly. For the nonlinear coefficient model it speci- 
fies the constant k 2 . The subscript IEQ varies from 1 to NEQ corresponding to 
each coupled equation. (See Table 3-4 for the order of the equations being solved.) 
Good values for a given application are usually determined by experience, but re- 
commended starting values are 1.0 for the constant coefficient model, 0.01 for the 
nonlinear model for flows without shocks, and 0.1 for the nonlinear model for 
flows with shocks. The default values are all 1.0. 

CAVS2I Second-order implicit artificial viscosity coefficient, z,, specified as CAVS2I(IEQ). 

The subscript IEQ varies from 1 to NEQ corresponding to each coupled equation. 
(See Table 3-4 for the order of the equations being solved.) Good values for a 
given application are usually determined by experience, but recommended starting 
values are 2.0 for the constant coefficient model, and 0.0 for the nonlinear model. 
The default values are all 2.0. 


Time Difference Centering Parameters 

The following parameters specify the type of time differencing scheme to be used. The generalized Beam 
and Warming (1978) time differencing formula is given by equation (3.1) of Volume 1. 


THC A 2-element array specifying the time differencing centering parameters 0 } and 0 2 

to be used for the continuity equation. The default values are 1.0 and 0.0. 

THX A 3-element array specifying the time differencing centering parameters 0i, 0 2 , and 

0 3 to be used for the x- momentum equation. The default values are 1.0, 0.0, and 
0.0. 

THY A 3-element array specifying the time differencing centering parameters 0\, 02, and 

0 3 to be used for the j-momentum equation. The default values are 1.0, 0.0, and 

0 . 0 . 

THZ A 3-element array specifying the time differencing centering parameters 0i , 02, and 

0 3 to be used for the z-momentum equation. The default values are 1.0, 0.0, and 

0 . 0 . 

THE A 3-element array specifying the time differencing centering parameters 8 U 0 2 , and 

0 3 to be used for the energy equation. The default values are 1.0, 0.0, and 0.0. 


THKE A 2-element array specifying the time differencing centering parameters 0i and 0 2 

to be used for the k-t equations. The default values are 1.0 and 0.0. 

The following table summarizes the time differencing schemes that may be used. The Euler implicit 
method is recommended for steady flows, and the 3-point backward implicit method is recommended for 
unsteady flows. 


0j 0 2 0 3 Method 


Accuracy 


l 

1/2 

1 


0 0 Euler implicit 

0 1/2 Trapezoidal implicit 

1/2 1 3-point backward implicit 


O(At) 

O(At) 2 

0(At) 2 


32 3.0 Input Description 


Proteus 3-D User's Guide 



3.1 >9 Namelist TIME 


Time Step Selection Parameters 

The following parameters determine the procedure used to set the time step size for the mean flow 
equations, and to change it as the solution proceeds. The various options for IDTAU are described in more 
detail in the description of subroutine TIMSTP in Section 4.2 of Volume 3, 

IDTMOD The time step size At is recomputed every IDTMOD'th step. The default value 
is 100000. 


IDTAU 


1 for a global (i.e., constant in space) time step At = (CFL)At c/7 , where At c/I is the 
minimum of the allowable time steps at each grid point based on the CFL cri- 
teria for explicit methods. 

2 for a global time step initially computed using the IDTAU = 1 option, but ad- 
justed as the solution proceeds based on AQ mo ,, the absolute value of the 
maximum change in the dependent variables. 8 For any of the dependent vari- 
ables, if AQ moi < CHG1, the CFL number is multiplied by DTF1. If 
AQ„„ > CHG2, the CFL number is divided by DTF2. If AQ„„ >0.15, the 
CFL number is cut in half. The CFL number will not be decreased below 
CFLMIN, or increased above CFLMAX. 

3 for a global time step At equal to the specified input DT. 

4 for a global time step initially equal to the specified input DT, but adjusted as 
the solution proceeds based, on AQ m „, the absolute value of the maximum 
change in the dependent variables. 8 For any of the dependent variables, if 
AQ™, < CHG1, At is multiplied by DTF1. If AQ m „ > CHG2, At is divided 
by DTF2. If AQ m „ >0.15, At is cut in half. At will not be decreased below 
DTMIN, or increased above DTMAX. 

5 for a local (i.e., varying in space) time step At = (CFL)At,/,, where A? c /i is the 
allowable time step at each grid point based on the CFL criteria for explicit 
methods. 

6 for a local time step initially computed using the IDTAU = 5 option, but ad- 
justed as the solution proceeds based on AQ mo ,, the absolute value of the 
maximum change in the dependent variables. 8 For any of the dependent vari- 
ables, if AQ™, < CHG1, the CFL number is multiplied by DTF1. If 
AQ m „ > CHG2, the CFL number is divided by UTF2. If AQ m „ >0.15, the 
CFL number is cut in half. The CFL number will not be decreased below 
CFLMIN, or increased above CFLMAX. 

7 for a global time step with cycling. At will be cycled repeatedly between 
DTMIN and DTMAX using a logarithmic progression over NDTCYC time 
steps. For some problems this option has been shown to dramatically speed 
convergence. However, the choice of DTMIN, DTMAX, and NDTCYC is 
critical, and no method has been developed that assures a good choice. Poor 
choices may even slow down convergence, so this option should be used with 
caution. 

8 for a local time step computed using the procedure of Knight and Choi (1989). 

With this option, At = max[AT 0 , (At c;7 ),]. where 

At 0 = (CFL) min[(AT c/ ,) { , (At cfl \, (At c/7 ) 5 ]. The parameters (At c// ) {> (At,/,),, and 
(At ,/j) c are the allowable time steps at each grid point based on the CFL criteria 
for explicit methods, computed separately for each computational coordinate 
direction. This formulation assumes that flow is generally in the £ direction. 

9 for a local time step computed as in the IDTAU = 8 option, but with a viscous 
correction added to the definitions of (At,/,) ,, (At ,/,),, and (At.,-),, similar to that 
used by Cooper (1987). 


8 In A Qmcx, the total energy Ej has been divided by Er, = p r RT r \(y r — l) + p r u?j2 so that it is the same order of 
magnitude as the other conservation variables. 
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If IDTAU — 7, ICHECK and IDTMOD are both automatically set equal to 1, and 
NITAVG is set equal to NDTCYC. In addition, if IDTAU = 7 and 
ICTEST = 1, ICTEST is changed to 2. If IDTAU = 2, 4, or 6, IDTMOD is au- 
tomatically set equal to ICHECK. The default value is 5. 

The above parameters IDTAU and IDTMOD apply to every case. Which of the remaining parameters 
are needed depends on the value of IDTAU, as specified in the following table. 

IDTAU Parameters Needed 

1 CFL 

2 CFL, CHG1, CHG2, DTF1, DTF2, CFLMIN, CFLMAX 

3 DT 

4 DT, CHG1, CHG2, DTF1, DTF2, DTMIN, DTMAX 

5 CFL 

6 CFL, CHG1, CHG2, DTF1, DTF2, CFLMIN, CFLMAX 

7 DTMIN, DTMAX, NDTCYC 

8 CFL 

9 CFL 

CFL An array, given as CFL(ITSEQ), specifying the ratio At/At c/v , where At is the ac- 

tual time step used in the implicit calculation and Ar c/I is the allowable time step 
based on the CFL criteria for explicit methods. The subscript ITSEQ is the se- 
quence number. For time steps 1 through NTIME(I), CFL(l) will be used. Then 
for steps NTIME(l) + 1 through NTIME(I) + NTIME(2), CFL(2) will be used, 
etc. 9 CFL is only used if IDTAU = 1, 2, 5, 6, 8, or 9. The default values are all 
1.0. 

DT An array, given as DT(ITSEQ), specifying the time step size At. The subscript 

ITSEQ is the sequence number. For time steps 1 through NTIME(l), DT(1) will 
be used. Then for steps NTIME(l) + I through NTIME(l) + NTIME(2), DT(2) 
will be used, etc. 10 DT is only used if IDTAU = 3 or 4. The default values are all 
0 . 01 . 

CHG1 Minimum change, in absolute value, that is allowed in any dependent variable 

before increasing At. CHG1 is only used if IDTAU = 2, 4, or 6. The default value 
is 0.04. 

CHG2 Maximum change, in absolute value, that is allowed in any dependent variable 

before decreasing At. CHG2 is only used if IDTAU = 2, 4, or 6. The default value 
is 0.06. 

DTF 1 Factor by which At is multiplied if the solution changes too slowly. DTF 1 is only 

used if IDTAU = 2, 4, or 6. The default value is 1.25. 

DTF2 Factor by which At is divided if the solution changes too quickly. DTF 2 is only 

used if IDTAU = 2, 4, or 6. The default value is 1.25. 

CFLMIN Minimum value that the CFL number is allowed to reach. CFLMIN is only used 
if IDTAU = 2 or 6. The default value is 0.5. 


9 Note that if IDTAU = 2 or 6, CFL(l) only sets At for the first time step, and that the time step sequencing option 
does not apply. 

>o Note that if IDTAU = 4, DT(1) only sets At for the first time step, and that the time step sequencing option does 
not apply. 
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CFLMAX Maximum value that the CFL number is allowed to reach. CFLMAX is only used 
if IDTAU = 2 or 6. The default value is 10.0. 


DTMIN 


DTMAX 


NDTCYC 


Minimum value that At is allowed to reach (IDTAU = 4), or the minimum At in 
the time step cycling procedure (IDTAU = 7.) The default value is 0.1. 

Maximum value that At is allowed to reach (IDTAU = 4), or the maximum At 
in the time step cycling procedure (IDTAU = 7.) The default value is 0.1. 

Number of time steps per time step cycle. NDTCYC is used only with 
IDTAU = 7. The default value is 2, which results in a constant At if 
DTMIN = DTMAX. 


Time Marching Limits 

These parameters determine the maximum number of time steps that will be taken. 

NTSEQ The number of time step sequences being used. The maximum value allowed is 

the value of the dimensioning parameter NTSEQP. If NTSEQ > 1, IDTAU must 
be equal to 1,3, or 5. (See Section 6.2.) The default value is 1. 

NTIME An array, given as NTIME(ITSEQ), specifying the maximum number of time 

steps to march. The subscript ITSEQ varies from I to NTSEQ, and allows a series 
of different time steps to be specified by the values of CFL or DT. 
NTIME(ITSEQ) specifies the number of time steps within sequence ITSEQ. If 
NTSEQ = 3, for example, the total number of time steps taken will be 
N tot ai = NTIME(l) -I- NTIME(2) + NTIME(3). The initial time level is level 1, 
and the final computed time level will be level N total + 1 . The default values are 
10, 9*0. 

TLIM When the amount of CPU time remaining for the job drops below TLIM seconds, 

the calculation is stopped. The default value is 20.0. 


Convergence Testing Parameters 

These parameters determine the convergence criteria to be used. 


ICHECK Results are checked for convergence every ICHECKth time level. The default 
value is 10. 


ICTEST 


1 to determine convergence based on the maximum change in absolute value of 
each of the conservation variables over a single time step, AQ mox . u 

2 to determine convergence based on the maximum change in absolute value of 
each of the conservation variables, averaged over the last NITAVG time steps, 
AQav*. 11 

3 to determine convergence based on Rl 2 , the L* norm of the residual for each 
equation. 

4 to determine convergence based on R^, the average absolute value of the resi- 
dual for each equation. 

5 to determine convergence based on R max , the maximum absolute value of the 
residual for each equation. 

Convergence is assumed when the maximum change or residual parameter is less 

than EPS. Note that the change in conservation variables over a time step is di- 
rectly related to the size of the time step. Small time steps naturally yield small 


11 The total energy Er is divided by Ej r = p r R T/(» — 1) + p r u}j2 before testing for convergence, so that it is the same 
order of magnitude as the other conservation variables. 
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changes in conservation variables. With ICTEST = 1 or 2, therefore, convergence 
may be indicated prematurely. 

If ICTEST = 2, ICHECK and IDTMOD are automatically set equal to 1. The 
default value is 3. 

EPS Level of convergence to be reached, specified as EPS(IVAR) where IVAR varies 

from 1 to NEQ, corresponding to each conservation variable or equation. The 
default values are all 0.001. 

NITAVG Number of time steps over which the maximum change in conservation variables 
is averaged to determine convergence. The maximum value allowed is the value 
of the dimensioning parameter NAMAX. (See Section 6.2.) NITAVG only ap- 
plies to the ICTEST = 2 option. The default value is 10. 


3.1.10 Namelist TURB 


Control Parameters 

The following parameters determine the type of turbulence model that will be used, and where it will 
be applied. These parameters apply to both the Baldwin-Lomax algebraic model and the Chien two- 
equation model. 

ITURB 0 for laminar flow. 

1 for turbulent flow, using the algebraic eddy viscosity model of Baldwin and 
Lomax (1978), as described in Section 9.0 of Volume 1. 

20 for turbulent flow, using the two-equation k-t model of Chien (1982), as de- 
scribed in Section 9.0 of Volume 1. 

The default value is 0. 


PRT If PRT > 0.0, it specifies the turbulent Prandtl number, which will be treated as 

constant. If PRT < 0.0, the turbulent Prandtl number will vary, and be computed 
using the empirical formula of Wassel and Catton (1973). The default value is 
0.91. 

LWALL1 A three-dimensional array, given as L3VALL 1(12, 13, IBOUND), specifying which 
points on the boundaries are on a solid wall. Here 12 = 1 to N2 and 13 = 1 to 
N3 corresponding to each grid point on the boundary, and IBOUND = I or 2, 
corresponding to the £ = 0 and £ = 1 boundaries, respectively. Valid values of 
LWALL1(I2, 13, IBOUND) are: 

0 if the 12,13 point on boundary IBOUND is not on a solid wall. 

1 if the 12,13 point on boundary IBOUND is on a solid wall. 

If LWALL1 is not specified, wall points on the £ boundaries are identified by 
looking for points with zero velocity. Therefore, LWALL1 need only be specified 
for cases with non-zero velocity at the wall (e.g., a moving wall, bleed, blowing, 
etc.). LWALL1 is not needed if the boundary condition for boundary IBOUND 
is set using the KBC1 (IBOUND) meta flag. The default values are all 0. 

LWALL2 A three-dimensional array, given as LWALL2(1 1,13, IBOUND), specifying which 
points on the rj boundaries are on a solid wall. Here 11 = 1 to N 1 and 13 = 1 to 
N3 corresponding to each grid point on the boundary, and IBOUND = 1 or 2, 
corresponding to the rj = 0 and = 1 boundaries, respectively. Valid values of 
LWALL2(I1, 13, IBOUND) are: 

0 if the 11,13 point on boundary IBOUND is not on a solid wall. 

1 if the 11,13 point on boundary IBOUND is on a solid wall. 
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If LWALL2 is not specified, wall points on the rj boundaries are identified by 
looking for points with zero velocity. Therefore, LWALL2 need only be specified 
for cases with non-zero velocity at the wall (e.g., a moving wall, bleed, blowing, 
etc.). LWALL2 is not needed if the boundary condition for boundary IBOUND 
is set using the KBC2(IBOUND) meta flag. The default values are all 0. 

LWALL3 A three-dimensional array, given as LWALL3(1 1,12, IBOUND), specifying which 
points on the £ boundaries are on a solid wall. Here 11 = 1 to NT and 12= 1 to 
N2 corresponding to each grid point on the boundary, and IBOUND = 1 or 2, 
corresponding to the £ = 0 and £ = I boundaries, respectively. Valid values of 
LWALL3(1 1,12, IBOUND) are: 

0 if the 11,12 point on boundary IBOUND is not on a solid wall. 

1 if the 11,12 point on boundary IBOUND is on a solid wall. 

If LWALL3 is not specified, wall points on the £ boundaries are identified by 
looking for points with zero velocity. Therefore, LWALL3 need only be specified 
for cases with non-zero velocity at the wall (e.g., a moving wall, bleed, blowing, 
etc.). LWALL3 is not needed if the boundary condition for boundary IBOUND 
is set using the KBC3(IBOUND) meta flag. The default values are all 0. 

Baldwin- Lomax Turbulence Model Parameters 

These parameters apply only to the Baldwin-Lomax algebraic eddy viscosity model. Note, however, 
that they are used wLen the Baldwin-Lomax model is used to generate initial turbulent viscosity values for 
the Chien k-z model (i.e., when IREST = 0, 1, or 3.) 

The following two options may be used to modify the Baldwin-Lomax model. 

INNER 1 to use the inner layer model of Baldwin and Lomax (1978). 

2 to use the inner layer model of Spalding (1961) and Kleinstein (1967). 

The default value is 1. 

ILDAMP 0 to use the normal Baldwin-Lomax mixing length formula in the inner region. 

1 to use the modified mixing length formula of Lau der and Priddin (1973) in the 
inner region of the Baldwin-Lomax model. 

The default value is 0. 

The following parameters are various constants used in the Baldwin-Lomax model. 

The Clauser constant K used in the Baldwin-Lomax outer region model. The de- 
fault value is 0.0168. 

The constant C cp used in the Baldwin-Lomax outer region model. The default 
value is 1.6. 

The constant C w k used in the formula for F wa ke in the Baldwin-Lomax outer region 
model. The default value is 0.25. 

The constant C^ub used in the formula for the Klebanoff intermittency factor F^ub 
in the Baldwin-Lomax outer region model. The default value is 0.3. 

The constant ( C K!eb ) min used in the formula for the Klebanoff intermittency factor 
F}ud> in the Baldwin-Lomax outer region model. The default value is normally 0.0. 
However, when the Baldwin-Lomax model is being used to generate initial turbu- 
lent viscosity values for the Chien k-t model (ITURB = 20 and IREST = 0, 1, or 
3), the default value is 0.1. 


CCLAU 

CCP 

CWK 

CKLEB 

CKMIN 
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APLUS The Van Driest damping constant A" used in the Baldwin- Lomax outer and inner 

region models. The default value is 26.0. 

CB The constant B used in the formula for the Klebanoff intermittency factor Fxieb in 

the Baldwin- Lomax outer region model, and in the Spalding- Kleinstein inner re- 
gion model. The default value is 5.5. 

CVK The Von K arm an mixing length constant k used in both the Baldwin- Lomax and 

Spalding- Kleinstein inner region models. The default value is 0.4. 

CNL The exponent n in the Launder-Priddin modified mixing length formula. The de- 

fault value is 1.7. 

Chien Turbulence Model Parameters 


These parameters apply only to the Chien k-z model. The following two parameters specify the time 
step size for the k-z equations. 

NTKE Number of k-z time iterations per mean flow time iteration. Note that this must 

be an integer equal to or greater than 1. Past experience indicates that this value 
should generally be less than 50. The default value is 1. 

TFACT Factor used in computing the k-z time step. Note that this can be any real positive 

number. The time step size for the k-z equations will be equal to 
At*., = TFACT(At), where At is the time step size for the mean flow equations. 
For equal k-z and mean flow time steps, use TFACT = 1/NTKE. The default 
value is 1.0. 


The following parameters are various constants used in the Chien k-z model. 


CMUR 

CONE 


Constant C Mr used to compute C„ in the turbulent viscosity formula for the k-z 
equations. The default value is 0.09. 

Constant Q used in the production term of the z equation. The default value is 
1.35. 


CTHREE 

CTWOR 

SIGE 

SIGK 


Constant C 3 used to compute C M in the turbulent viscosity formula for the k-z 
equations. The default value is 0.0115. 

Constant C 2r used to compute C 2 in the dissipation term of the z equation. The 
default value is 1.8. 

The constant a, used in the diffusion term of the z equation. The default value is 
1.3. 

The constant o k used in the diffusion term of the k equation. The default value is 

1.0. 


3.1.11 Namelist IC 

This namelist is used in subroutine INIT to read in parameters needed in setting up the initial conditions 
for the mean flow. The version of INIT built into Proteus specifies uniform flow with constant properties 
everywhere in the flow field. In general, however, the user will supply a version of INIT tailored to the 
problem being solved. This namelist, then, may be modified by the user to read in parameters different from 
those listed here. 


P0 Initial static pressure /%. The default value is 1.0. 

TO Initial static temperature T 0 . The default value is 1.0. 
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uo 

Initial ^-direction velocity Uq . 

The default value is 0.0. 

VO 

Initial ^-direction velocity v 0 . 

The default value is 0.0. 

wo 

Initial z-direction velocity w 0 . 

The default value is 0.0. 

COORDINATE SYSTEM FILE 



The type of computational coordinate system to be used is controlled by the input parameter NGEOM 
in namelist GMTRY. For NGEOM = 10, the coordinate system is read from a pre-stored file. This file 
may be created by any body-fitted coordinate system generator available to the user. The coordinates may 
be nonorthogonal. 

The metric coefficients and Jacobian describing the nonorthogonal grid transformation are computed 
internally by Proteus . This calculation involves numerically computing first derivatives of the user- specified 
coordinates. Since Proteus solves the Navier-Stokes equations in fully conservative form, the metric coef- 
ficients themselves are factors in terms whose first and second derivatives are also computed numerically. 
In effect, then, third derivatives of the user-specified coordinates are used in the solution. Care should 
therefore be taken in ensuring that these coordinates are smooth. No coordinate smoothing is done by 
Proteus itself. 

The Cartesian (;c j>,z) coordinates describing the computational coordinate system are read from an un- 
formatted file as follows: 

read (ngrid) ngl,ng2,ng3 

read (ngrid) ( ( (xc( jl , 52, j3) , jl = l >ngl ) , 52-\ >ng2) , j3=l ,ng3) , 

$ ( ( (yc( jl , 52, j3) , jl = l ,ngl) , j2=l ,ng2) , J3=l ,ng3) , 

$ C((zc(jl, 32 ,j 3 ),.jl = l,ngl), j2-l ,ng2) , j3-l ,ng3) 

The parameters read from the file are defined as follows: 

NG1 Number of points in the { direction. The maximum value allowed is the value of 

the dimensioning parameter NIP. (See Section 6.2.) 

NG2 Number of points in the rj direction. The maximum value allowed is the value of 

the dimensioning parameter N2P. (See Section 6.2.) 

NG3 Number of points in the { direction. The maximum value allowed is the value of 

the dimensioning parameter N3P. (See Section 6.2.) 

XC Cartesian x-coordinate. 

YC Cartesian j-coordinate. 

ZC Cartesian z-coordinate. 

Note that the number of points NG1, NG2, and NG3 used to specify the computational coordinate 
system need not be the same as the number of points Nl, N2, and N3 used in the computational mesh . 
The coordinates of the points in the computational mesh, which is the mesh used in the Proteus solution, 
will be found by interpolation among the points in the computational coordinate system. 
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TABLE 3-1. - NORMALIZING CONDITIONS 


Variable 

Normalizing Value 

Length 

Ln = U 

Velocity 


Temperature 

T n =T, 

Density 

Pn = Pr 

Viscosity 


Thermal conductivity 

k n = k r 

Pressure 

> 

II 

Energy per unit volume 

e n - p r u} 

Gas constant 

R* - u?jT r 

Specific heat 

Cpn = tfIT, 

Enthalpy 

K, = U} 

Time 

t„ - LrlUr 

Turbulent kinetic energy 

k„ = u? 

Turbulent dissipation rate 

£n = prl^jPr 


TABLE 3-2. - REFERENCE CONDITIONS 


Variable 

Reference Value 

Length 

Lr 

Velocity 


Temperature 

T r 

Density 

Pr 

Viscosity 

Pr 

Thermal conductivity 

K_ 

Pressure 

Pr = p,R Trlgc 

Energy per unit volume 

e r = p r u} 

Enthalpy 

U? 

Specific heat 

U?ITr 

Time 

L/^r 

Turbulent kinetic energy 

u} 

Turbulent dissipation rate 

PrUflPr 
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TABLE 3-3. - OUTPUT VARIABLES 


IVOUT 

VARIABLE 

DEFINITION 

Velocities 

i 

x-velocity 

u 

2 

y-velocity 

V 

3 

z-velocity 

w 

4 

Mach number 

„ 11 
M= a 

5 

Speed of sound 

a = JyRT 

6 

Contravariant velocity 
normal to £ surface 

a , + ui x + v£ y + +f y + g)' 1 12 

7 

Contravariant velocity 
norma] to r\ surface 

{f] l + uri x + vr] y + wr} 2 )l(rj x + rjj + »J^)’ /2 

K, + + <y + wC z )(t; 2 x + C y + si) 1/2 

8 

Contravariant velocity 
normal to f surface 


9 

Total velocity magnitude 

| V | = (u 2 -F v 2 + w 2 y 

10 

x-momentum 

pu 

11 

y-momentum 

pv 

12 

z-momentum 

pw 

13 

£ -velocity 

= (x^u + y^v + z ? w)/(xj + y\ + Zj) 1/2 

14 

^-velocity 

U, = {x n u + + z^)]{x\ + y 2 + z 2 )' 12 

15 

£■ -velocity 

F c = (x^u + y^v + z^w)l(xf + y 2 + z 2 )' 12 

16 

Flow angle a„, deg. 

tan 

17 

Flow angle a W} deg. 

tan_ 1 17 

Densities 

20 

Static density 

P 

/ i \ */(y- 1) 

Pr= + 2 ) 

21 

Total density 
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IVOUT 

VARIABLE 

DEFINITION 

Pressures 


Static pressure 
Total pressure 

Static pressure coefficient 

Total pressure coefficient 
Pitot pressure 


p T = p(l+^r—M 2 


P -Pr 

C p 2 n 

Pr u r ! 2 &c 

PT ~ PT r 

CpT ~ Pr u}l2g c 
Pp=PT if a/ < 1 


y Hr- i) 


p p~p{ r r~ M2 ) 


y l(y- i) 


y 4- 1 y 4- 1 


— l/(y— 1) 


if M> 1 


Dynamic pressure 


Static temperature 
Total temperature 


I r 2 , 2 , . ,2\ Pr u r 

- P (u +V +w) — 


Temperatures 

T 

t t =t(\ + 


Energies 


y - 1 t ,2 

2 M 


50 

Total energy per unit volume 

E T 

51 

Total energy 

E r 

p 

52 

Internal energy 

e i = c v T 

53 

Kinetic energy 

e k = \{u 2 + v 2 + w 2 ) 

Enthalpies 

60 

Static enthalpy 

II 

61 

Total enthalpy 

h T — c p T T 
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IVOUT 

VARIABLE 

DEFINITION 

Vorticities 

70 

je-vorticity 

0 dw dv 

x ~ dy dz 

71 

j/-vorticity 

„ du dw 

y - dz dx 

72 

z-vorticity 

n dv du 

z ~ dx dy 

73 

Total vorticity magnitude 

|n| =(n 2 x + n 2 y + n ]) li2 

Entropies 

80 

Entropy 


Temperature- Dependent Parameters 

90 

Laminar viscosity coefficient 

=t 

I 

ii 

£ 

91 

Laminar second coefficient of 

2 M 


viscosity 

X >~ 3 

92 

Laminar thermal conductivity 
coefficient 

1 

II 

93 

Specific heat at constant 
pressure 

C P 

94 

Specific heat at constant vol- 

Cy 


ume 


95 

Ratio of specific heats 

II 
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IVOUT 


VARIABLE 


DEFINITION 


Turbulence Parameters 


100 

Turbulent viscosity coeffi- 
cient 

Pi 

101 

Turbulent second coefficient 
of viscosity' 

, 2 Pt 

kt ~ T 

102 

Turbulent thermal 
conductivity coefficient 

i 

' Pr t K 

103 

Effective viscosity coefficient 

M 

104 

Effective second coefficient 
of viscosity 

A 

105 

Effective thermal 
conductivity coefficient 

k 

106 

Turbulent kinetic energy 

k 

107 

Turbulent dissipation rate 

\ P^y* n . 

T + - fJL w Re ' 

108 

Inner region coordinate 


109 

Inner region velocity 

Gradients 


120 

Shear stress 


121 

Shear stress 

Tyy 

122 

Shear stress 


123 

Shear stress 

Txy 

124 

Shear stress 

Txz 

125 

Shear stress 

A* 

126 

Heat flux 

<h 

127 

Heat flux 

q> 

128 

Heat flux 
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IVOUT 

VARIABLE 

DEFINITION 

Coordinates 

200 

Cartesian ^-coordinate 

X 

201 

Cartesian ^-coordinate 

y 

202 

Cartesian z-coordinate 

Z 

203 

Local Ax 

l-T’l [(n/z - nly) Ac + - ZyQAri 




204 

Local A y 

\J~ 1 1 lto£ x ~ + {Uz ~ «/*) M 




205 

Local Az 

| r 1 1 [inly - «!/*) + (*A - Uy)^n 



+ (Zx*l y - iyVx) A i] 

206 

Local cell Reynolds number 
Re c 

p|u|[(Ax) 2 + (A^) 2 + (Az) 2 ] 1/2 
Re r ^ 

207 

Local x-direction cell 
Reynolds number ( Re c ) x 

p | «| Ax 

208 

Local ^-direction cell 
Reynolds number ( Re c ) y 

■ n pMAj; 
Re r p 

209 

Local z-direction cell 
Reynolds number ( Re c ) z 

p | w | Az 

^ e r p 

Metric Parameters 

210 

Inverse of the grid transfor- 
mation Jacobian 

r x 

211 

Metric coefficient 

it 

212 

Metric coefficient 

Zx 

213 

Metric coefficient 

iy 

214 

Metric coefficient 

iz 

215 

Metric coefficient 

Vt 

216 

Metric coefficient 

Vx 

217 

Metric coefficient 

*iy 

218 

Metric coefficient 

n 2 

219 

Metric coefficient 


220 

Metric coefficient 

ix 

221 

Metric coefficient 

iy 

222 

Metric coefficient 

iz 


Proteus 3-D User's Guide 


3.0 Input Description 45 






IVOUT VARIABLE 


DEFINITION 


Time step size 
Time 

Local CFL number 


Local ^-direction CFL num- 
ber 


Local ^-direction CFL num- 
ber 


Local C -direction CFL num- 
ber 


\ y h+ Ur 1 x + Vr ly + Wf lz I 

|C,+ «C x + vC y + wC 2 l 
+ AC 

f Cjc fjjc Cx 

+ a _ XT + ^r L+ ^F 
\ AC An AC 



+ -^ + 


AC 

A n + 

AC 

** 


k 

AC 

+ An + 

AC 


P 


+ U£ X + V£y + W£ z | 


(Cx 2 + ^/ + 0 1/2 


l»7r + “>/x + w /y + w »fel 
A>7 

(»?x 2 + »?/ + v z 2 y 12 
+ a ^ 

|C f + «Cx + V ^ + K z l 
AC 

(Cx 2 + C y 2 + C 2 2 ) 1/2 
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TABLE 3-4. - EQUATIONS SOLVED 


IHSTAG 

NEQ 

Order of Equations 

Order of 

Dependent Variables 

1 

4 

Continuity, x-momentum, 
^-momentum, z-momentum 

p, pu, pv, pw 

0 

5 

Continuity, x-momentum, 
^-momentum, z-momentum, energy 

p, pu, pv, pw, Et 
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TABLE 3-5. - BOUNDARY TYPES 


KBC 

VALUE 3 

BOUNDARY TYPE 

JBC 

VALUES SET 

EQUATIONS 

± 1 

No-slip adiabatic wall 

11, 21, 31, ±43, 
± 53 

u = v = w = 0, 8 pjdn — dTjdn = 0 

±2 

No -slip wall, specified tem- 
perature 

11, 21, 31, ±43, 
50 

u = v = w = 0, dpjdn — 0, AT = 0 

±3 

Inviscid wall 

±43, ±53, 71, 
m, 77* 

dpjdn = dTjdn = 0, V„ = 0, 
d 2 V tl jd<t> 2 = 0, d 2 Vald<t> 2 = 0 

10 

Subsonic inflow, linear ex- 
trapolation 

14, 24, 34, 46, 56 

d 2 ujd<t> 2 = d 2 vjd<f> 2 = d 2 wld<j> 2 = 0, 
Apr = A7> = 0 

± 11 

Subsonic inflow, zero gradi- 
ent 

+ 12, + 22, + 32, 
46, 56 

duld(f> = dvjd4> — dwjd<f> = 0, 
Ap r = ATj = 0 

20 

Subsonic outflow, linear ex- 
trapolation 

14, 24, 34, 40, 54 

d 2 ujd4> 2 = d 2 vjd(f> 2 = d 2 wjd<j> 2 = 0, 
A p = 0, d 2 Tjd4> 2 = 0 

±21 

Subsonic outflow, zero gradi- 
ent 

+ 12, + 22, + 32, 
40, ±52 

dujd<t> = dvjd(j> = dwjd<f> = 0, 
Ap = 0, dTjd<f> = 0 

30 

Supersonic inflow 

10, 20, 30, 40, 50 

Am = Av = Aw = Ap = AT = 0 

40 

Supersonic outflow, linear 
extrapolation 

14, 24, 34, 44, 54 

d 2 ujd<f> 2 = d 2 vjd<j> 2 = d 2 wjd<p 2 = 0, 
d 2 pjd<j> 2 = d 2 Tjd<j> 2 = 0 

±41 

Supersonic outflow, zero gra- 
dient 

4- 12, + 22, + 32, 
±42, ±52 

dujdcji = dvjd<j> = dwjd<f> = 0, 
dpjd<f> = dTjd<f> = 0 

±50 

Symmetry 

i 43, + 53, 71 , 
m, n c 

dpjdn = dTjdn = 0, V n = 0, 
dV.Jdj, = 0, dV a jd4> = 0 

1 60 

Spatially periodic 


Ql = Qni, Ql — QjV2j ° r Ql = Q>’3 


a Use the *+* sign for 2-point one-sided differencing of first derivatives, and the * sign for 3-point differencing of 
first derivatives. 

b The JBC values used are 84 and 89 for a \ boundary, 79 and 89 for an rj boundary, and 79 and 84 for a £ boundary. 

c The JBC values used are ± 83 and ± 88 for a £ boundary, ± 78 and ± 88 for an rj boundary, and ± 78 and ± 83 
for a £ boundary. 
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TABLE 3-6. - BOUNDARY CONDITION TYPES 


JBC OR IBC 
VALUE 1 

EQUATION 

DESCRIPTION 

Conservation Variable Boundary Conditions 

0 

AQ = 0 

No change from initial conditions. 

1 

Q =/ 

Specified conservation variable. 

+ 2 

dQld<f> =/ 

Specified coordinate direction gradient. 

±3 

dQldn =/ 

Specified normal direction gradient 5 . 

4 

d 2 Qld<f> 2 = 0 

Linear extrapolation. 

5 


Not used. 

6 


Not used. 

7 


Not used. 

8 


Not used. 

9 


Not used. 

x-Velocity Boundary Conditions 

10 

Aw = 0 

No change from initial conditions. 

11 

u=f 

Specified x- velocity. 

+ 12 

duid<t> =/ 

Specified coordinate direction gradient. 

+ 13 

dujdn =f 

Specified normal direction gradient 5 . 

14 

d 2 ujd<f> 2 = 0 

Linear extrapolation. 

15 


Not used. 

16 


Not used. 

17 


Not used. 

18 


Not used. 

19 


Not used. 

^-Velocity Boundary Conditions 

20 

Av = 0 

No change from initial conditions. 

21 

v=/ 

Specified ^-velocity. 

+ 22 

dvjd4> =/ 

Specified coordinate direction gradient. 

+ 23 

dv/dn —f 

Specified normal direction gradient 5 . 

24 

d 2 vld<j> 2 = 0 

Linear extrapolation. 

25 


Not used. 

26 


Not used. 

27 


Not used. 

28 


Not used. 

29 

tan" ] (vju) = f 

Specified flow angle in degrees. 

z-Velocity Boundary Conditions 

30 

> 

II 

o 

No change from initial conditions. 

31 

w=/ 

Specified z- velocity. 

+ 32 

dwld<f> =/ 

Specified coordinate direction gradient. 

+ 33 

dwjdn = f 

Specified normal direction gradient 5 . 

34 

d 2 wjdc{> 2 = 0 

Linear extrapolation. 

35 


Not used. 

36 


Not used. 

37 


Not used. 

38 


Not used. 

39 

tan _1 (w/w) = / 

Specified flow angle in degrees. 
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4-1 +1 


JBC OR IBC 


EQUATION 


DESCRIPTION 


VALUE* 


Pressure Boundary Conditions 


40 

o 

II 

a 

No change from initial conditions. 

41 

p=f 

Specified static pressure. 

+ 42 

dpjdt =/ 

Specified coordinate direction gradient. 

+ 43 

dpjdn = / 

Specified normal direction gradient 11 . 

44 

d 2 pjd4> 2 = 0 

Linear extrapolation. 

45 


Not used. 

46 

A p T - 0 

No change from initial conditions. 

47 

Pr-f 

Specified total pressure. 

48 


Not used. 

49 


Not used. 


Temperature Boundary Conditions 


50 

> 

II 

O 

No change from initial conditions. 

51 

T—f 

Specified static temperature. 

52 

dTjd4> =f 

Specified coordinate direction gradient 

53 

dTjdn = f 

Specified normal direction gradient 5 . 

54 

PTIdt 2 = 0 

Linear extrapolation. 

55 

Not used. 

56 

O 

II 

£ 

< 

No change from initial conditions. 

57 

II 

Specified total temperature. 

58 


Not used. 

59 


Not used. 


Density Boundary Conditions 


60 

Ap — 0 

No change from initial conditions. 

61 

p=f 

Specified static density. 

+ 62 

dpjd<f> =/ 

Specified coordinate direction gradient. 

+ 63 

dpjdn = f 

Specified normal direction gradient 5 . 

64 

d 2 pjd4> 2 = 0 

Linear extrapolation. 

65 


Not used. 

66 


Not used. 

67 


Not used. 

68 


Not used. 

69 


Not used. 
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JBC OR IBC 
VALUE 1 

EQUATION 

DESCRIPTION 

Normal 0 and Coordinate Direction Velocity Boundary Conditions 

70 


Not used. 

71 

II 

Specified normal velocity. 

+ 72 

dv„m =/ 

Specified coordinate direction gradient. 

+ 73 

dV„/dn =/ 

Specified normal direction gradient 6 . 

74 

d 2 V„jd<t> 2 = 0 

Linear extrapolation. 

75 


Not used. 

76 

V t =f 

Specified f -velocity. 

+ 77 

dv,id<j> =f 

Specified coordinate direction gradient. 

+ 78 

dV t ldn—f 

Specified normal direction gradient 6 . 

79 

dWJdtf = 0 

Linear extrapolation. 

80 


Not used. 

81 

v*-f 

Specified ^-velocity. 

+ 82 

dVJd* =f 

Specified coordinate direction gradient. 

+ 83 

dVJdn-f 

Specified normal direction gradient 6 . 

84 

d 2 VJd4> 2 = 0 

Linear extrapolation. 

85 


Not used. 

86 

V(=f 

Specified £ -velocity. 

+ 87 

dvjd4> =/ 

Specified coordinate direction gradient. 

+ 88 

dVrldn=f 

Specified normal direction gradient 6 . 

89 

d 2 VJd<t> 2 = 0 

Linear extrapolation. 

User-Supplied Boundary Conditions 

90 

> 

II 

O 

No change from initial conditions. 

91 

F—f 

Specified function. 

+ 92 

dFjd<j} =/ 

Specified coordinate direction gradient. 

+ 93 

dFjdn =f 

Specified normal direction gradient 6 . 

94 

d 2 Fld<f> 2 = 0 

Linear extrapolation. 

95 


Not used. 

96 


Not used. 

97 


Not used. 

98 


Not used. 

99 


Not used. 


a Use the ' r -h' r sign for 2-point one-sided differencing, and the ' sign for 3-point one-sided differencing. 
b Normal derivatives are positive in the direction of increasing { or r\. 
c Normal velocity is positive in the direction of increasing { or rj. 


Proteus 3-D User's Guide 


3.0 Input Description 51 










TABLE 3-7. - BOUNDARY CONDITION TYPES FOR k-t EQUATIONS 


JBCT OR IBCT 
VALUE a 

EQUATION 

DESCRIPTION 

Turbulent Kinetic Energy Boundary Conditions 

0 

> 

II 

o 

No change from initial conditions. 

1 

k —f 

Specified turbulent kinetic energy. 

+ 2 

dkjd4> =f 

Specified coordinate direction gradient. 

3 


Not used. 

4 

d 2 kld<j> 2 = 0 

Linear extrapolation. 

5 


Not used. 

6 


Not used. 

7 


Not used. 

8 


Not used. 

9 


Not used. 

Turbulent Dissipation Rate Boundary Conditions 

10 

As = 0 

No change from initial conditions. 

11 

£=/ 

Specified turbulent dissipation rate. 

+ 12 

dtjd<j> =/ 

Specified coordinate direction gradient. 

13 


Not used. 

14 

dhjd4> 2 = 0 

Linear extrapolation. 

15 


Not used. 

16 


Not used. 

17 


Not used. 

18 


Not used. 

19 


Not used. 


* Use the '+' sign for 2-point one-sided differencing, and the ' sign for 3-point one-sided differencing. 
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4.0 OUTPUT DESCRIPTION 


Several output files may be created during a Proteus run. The standard output is a formatted file written 
to Fortran unit NOUT that is intended for printing. Additional unformatted files may be written for use 
as input by various post- processing programs. Unformatted restart files may also be written for use as input 
for a subsequent Proteus run. 

4.1 STANDARD OUTPUT 

The standard Proteus output is a formatted file written to Fortran unit NOUT, and is intended for 
printing. Actual examples of typical standard output files are presented in Section 9.0. Unless specified 
otherwise, all of the output parameters in the standard output are nondimensional, with the appropriate 
reference condition from Table 3-2 as the nondimensionalizing factor. 

4.1.1 Title Page and Namelists 

The standard Proteus output begins with a title page, 12 which identifies the version of Proteus being run 
and lists the user- specified title for the run. This is followed by a printout of the contents of the input 
namelists RSTRT, 10, GMTRY, FLOW, BC, NUM, TIME, and TURB. Note that, for variables not 
specified by the user in the input namelists, the values in this printout will be the default values. 

4.1.2 Normalizing and Reference Conditions 

The dimensional values for the normalizing and reference conditions are printed on the next page, with 
the appropriate units as set by the input parameter IUNITS. The normalizing conditions are the parameters 
used to nondimensionalize the governing equations. The reference conditions are used during input and 
output for nondimensionalization of various parameters and for specifying various flow conditions. The 
distinction between normalizing and reference conditions is described in greater detail in Section 3.1.1. They 
are listed in Tables 3-1 and 3-2. 

After the printout of the normalizing and reference parameters comes anything written to unit NOUT 
by the user-supplied subroutine INIT. For the default version of INIT supplied with Proteus, this consists 
only of the contents of namelist IC. 

4.1.3 Boundary Conditions 


The next page is a printout of the boundary conditions being used. The boundary condition parameters 
JBC and GBC are printed for the and £ boundaries. If the Chien two -equation turbulence model is 
being used, the boundary condition parameters for the k-i equations, JBCT and GBCT, are printed next. 
If time-dependent boundary conditions are being used, this is followed by a listing of the input tables of 
GTBC vs. NTBCA. 

4.1.4 Computed Flow Field 

The bulk of the standard Proteus output consists of printout of the computed flow field. The input array 
IVOUT determines which variables are printed, as described in Section 3.1.4. The variables currently 
available for printing are listed and defined in Table 3-3. The printout for each variable at a given tune le\el 
will begin on a separate page. The header for each variable will include the time level n , and, for global time 
steps (IDTAU =1-4, 7), the time t and time increment At in seconds. The flow field results are printed 


12 In this discussion, when 'pages' of output are referred to it is assumed that the file is printed with Fortran carriage 
control in effect 
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in blocks, with each block corresponding to a C location. Within each block, each column corresponds to 
a £ location, and each row to an t] location. The columns and rows are numbered with the l and rj indices. 

Flow field results are printed at time levels and grid points specified by the user through parameters in 
namelist 10. Since this printout can be very lengthy, the user is encouraged to minimize the amount of 
printed output by making judicious use of these parameters. Usually, the computed results can be examined 
most efficiently using post-processing graphics routines like CONTOUR or PL0T3D. (See Section 4.2). 

After the flow field printout, if the run ends normally a message is printed indicating whether or not the 
calculation converged. 

4.U5 Boundary Parameters 


After each flow field printout, various parameters may be printed along the boundaries. The input ar- 
rays IW0UT1, IW0UT2, and IW0UT3 determine whether or not these parameters are printed, as de- 
scribed in Section 3.1.4. The parameters printed are defined below. Note that some of these are meaningful 
only if the boundary is a solid wall. 


X, Y,Z 

P 

CF 


TAUW 


T 

QW 


H 


ST 


Cartesian coordinates x, y, and z. 
Static pressure p . 

Skin friction coefficient C/, defined as 


c /= 


_ &V { 

1 2 


_ 2 _ 

Re r ^ dn 


where the overbar denotes a dimensional quantity. In this equation dV t jdn 
represents the normal derivative of the tangential velocity. 

Shear stress t„, defined as 


t* is thus nondimensionalized by 
Static temperature T . 

Heat flux q W7 defined as 



In this equation dTfdn represents the normal derivative of the temperature. q w 
is thus nondimensionalized by k r T r jLr. 

Heat transfer coefficient h , defined as 




T — 1 


h is thus nondimensionalized by k r \Lr. 
Stanton number St } defined as 


St = 


h 

PPfip 


h 1 
c p Re r Pr r 


54 4.0 Output Description 


Proteus 3-D User's Guide 



The boundary parameters are printed at the same time levels as the flow field printout, and at the points 
on the boundary' specified by the input parameters IPRT1, IPRT2, and IPRT3, or IPRT1A, IPRT2A, and 
IPRT3A. 

For this printout, normal derivatives are computed using a normal vector n directed into the flow field. 
This means, for example, that the skin friction C/ will be positive for attached flow- , even on the upper 
boundary. 

4.1.6 Convergence History 

In evaluating the results of a steady Proteus calculation, it's important to consider the level of conver- 
gence. This may be done by examining one of the forms of the residual for each equation. The residual 
is simply the number resulting from evaluating the steady form of the equation at a specific grid point and 
time (or iteration) level. Ideally, the residuals would all approach zero at convergence. In practice, however, 
for Teal" problems they often drop to a certain level and then level off. Continuing the calculation beyond 
this point w'ill not improve the results. 

A decrease in the L 2 norm of the residual of three orders of magnitude is sometimes considered sufficient. 
Convergence, however, is in the eye of the beholder. The amount of decrease in the residual necessary for 
convergence will vary from problem to problem. For some problems, it may even be more appropriate to 
measure convergence by some flow-related parameter, such as the lift coefficient for an airfoil. Deteimining 
when a solution is sufficiently converged is, in some respects, a skill best acquired through experience. 

At the end of a Proteus calculation, if first-order time differencing and steady boundary conditions were 
used, a summary of the convergence history is printed for each governing equation. 13 The parameters in this 
printout are defined as follows: 

LEVEL 
CHGMAX 


CHGAVG 


RESL2 


RESAVG 

RESMAX 

LRMAX 


Time level n. 

Maximum change in absolute value of the dependent variables from time level 

72 — 1 tO 77. 14 

' AQ majc = max |AQ" y ‘*| 

Maximum change in absolute value of the dependent variables, averaged over 
the last NITAVG time steps. 14 

AQav »° = NITAVG 

m= n -NITAVG 

The Li norm of the residual at time level n. 

The average absolute value of the residual at time level 77, R^. 

The maximum absolute value of the residual at time level n, R m „ 

The grid indices (/,/, k) corresponding to the location of R mox . 



13 Second-order time differencing should be used only for unsteady problems, for which "convergence' has no 
meaning. It should also be noted that the computation of the residuals in the code is correct only for first-order 
time differencing. 

14 For the energy equation, the change in Ej is divided by Ej r = p r R T r l(y r — 1) + Pru}j2 y so that it is the same order 
of magnitude as the other conservation variables. 
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In computing the residuals, the summations, maximums, and averages are over all interior grid points, plus 
points on spatially periodic boundaries. 

To avoid undesirably long tables, the convergence parameters are printed at an interval that limits the 
printout to NHMAX time levels. N'HMAX can be specified by the user in namelist 10. However, the 
residuals are always printed at the first two time levels. This is done because the residuals at tune level 1 
(the initial condition level) may not be truly representative of the degree of convergence. For instance, if 
the initial conditions are zero velocity and constant pressure and temperature at every interior point, the 
computed residuals will be exactly zero. When the time marching procedure begins, however, the flow field 
will start developing in response to the boundary conditions, and the residuals will reach a maximum in the 
first few time steps. Note that, in the printout, CHGMAX will be zero until time level n = ICHECK. 
CHGAVG will only be computed when ICTEST = 2, and will be zero until time level n = NITAVG. 

As noted in Section 8. 1 of Volume 1 , adding artificial viscosity changes the original governing partial 
differential equations. For cases run with artificial viscosity, therefore, the residuals are printed both with 
and without the artificial viscosity terms included. This may provide some estimate of the overall error in 
the solution introduced by the artificial viscosity. Convergence is determined by the residuals with the ar- 
tificial viscosity terms included. 

4.1.6 Additional Output 

In addition to the output discussed above, various types of additional printout can be generated by the 
IDEBUG options, as discussed in Section 3.1.4. Various diagnostics may also appear in the standard out- 
put file. These are discussed in greater detail in Section 7.0. 

4.2 PLOT FILES 

The amount of flow field data generated by a Navier- Stokes code is normally much too large to effi- 
ciently comprehend by examining printed output. The computed results are therefore generally examined 
graphically using various post-processing plotting routines. These plotting routines require as input a file 
or files, generally called plot files, that are written by the flow solver and contain the coordinates and com- 
puted flow field data. 

Various types of unformatted plot files may be written by Proteus, as controlled by the input parameter 
T FLOT discussed in Section 3.1.4. These files are designed for use by either the r ONTOUR or PLOT3D 
plotting programs. 15 

CONTOUR is a three-dimensional plotting program developed at NASA Lewis using internal Lewis- 
developed graphics routines. It currently can be used only at NASA Lewis on the Amdahl 5860 computer 
using the VM operating system, or from the Scientific VAX Cluster using the VMS operating system. 
Originally designed for use with three-dimensional Parabolized Navier- Stokes (PNS) codes, it can be used 
to generate various types of contour and velocity vector plots in computational planes. 

PLOT3D (Walatka, Buning, Pierce, and Elson, 1990) is a sophisticated three-dimensional plotting 
program specifically designed for displaying results of computational fluid dynamics analyses. It is widely 
used in government . industry, and universities for interactive visualization of complex flow field data gen- 
erated by CFD analyses. The computational grid is stored in one file, called an XYZ file, and the computed 
flow field is stored in another file, called a Q file. There are several options within PLOT3D concerning the 
format of these files. At NASA Lewis, PLOT3D is available on the Silicon Graphics IRIS workstations 
and on the Scientific VAX Cluster. 

It should be noted that, in Fortran, unformatted files are not generally transportable from computer to 
computer. If the plot files written by Proteus are to be used on some other computer (e.g., a graphics 
workstation), a separate conversion program will normally be required. 


•s If only the last computed time level is of interest, the restart files may also be used for plotting with PLOT3D. See 
Section 4.4. 
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4.2.1 CONTOUR Plot File (IPLOT = 1) 

With the IPLOT = 1 option in Proteus, a plot file is generated for use by CONTOUR using the fol- 
lowing Fortran statements: 

write Cnplot) title 

write Cnplot) machr , rer , lr , ur ,pr , tr , rhor , rg , gamr 
do 20 il = l,nl 

write C nplot) level , n2 , n3 , isym, system 
do 10 i2 = l,n2 _ . 

write Cnplot) ( (f (ivar /il > i2, 13) , ivar- 1 , 14) , i3-l >n3) 

10 continue 

20 continue 

All of the above WRITE statements are executed for each time level written into the file. The plot file thus 
consists of multiple sets of data, each containing the computed results at a single time level. Note that the 
value of the time t is not written into the file. Note also that Nl, the number of grid points in the £ di- 
rection, is not written into the file. 

Unless specified otherwise, all of the parameters written into the CONTOUR plot file are nondimen 
sional, with the appropriate reference condition as the nondimensionalizing factor. The parameters are de- 
fined as follows: 


TITLE 

MACHR 

RER 

LR 

UR 

PR 

TR 

RHOR 

RG 

GAMR 

LEVEL 

N2 

N3 

ISYM 

SYSTEM 


A descriptive title for the problem. 

Reference Mach number, M r = u r [Jy,RT r . This is a real variable. 
Reference Reynolds number, Re r = p r u r Lrjn r . 

Reference length L r in feet (meters). This is a real variable. 
Reference velocity u, in ft /sec (m/sec). 

Reference static pressure p r = p r RT r lg c in lb f /ft 2 (N/m 2 ). 

Reference temperature T r in °R (K). 

Reference density p r in lb m /ft 3 (kg/m 3 ). 

Gas constant R in ft 2 /sec 2 -°R (m 2 /sec 2 -K). 

Reference ratio of specific heats, y r = c^/<v 

Time level n. 

Number of grid points N 2 m the r\ direction. 

Number of grid points Ni in the £ direction. 

A symmetry parameter used in CONTOUR, set equal to 1. 

A coordinate system parameter used in CONTOUR, set equal to 0. 


F(l,„) 

Cartesian x coordinate. 

F(2,„) 

Cartesian y coordinate. 

F(3,„) 

Cartesian 2 coordinate. 

F(4,„) 

Velocity in the £ direction, V K . 

F(5,„) 

Velocity in the r\ direction, V r 

F(6,„) 

Velocity in the C direction, V c . 

F(7,„) 

Static pressure plp r . 

F(8,„) 

Static temperature T. 

F(9,„) 

Mach number M = Ju 2 + v 2 H- w 2 jJyRT 
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F(10,„) 

x-velocity u . 

F(1 1,„) 

^-velocity v. 

F(12,„) 

z- velocity w. 

F(13,„) 

Static density p. 

F(14,„) 

Total energy per unit volume E T . 

4.2.2 PLOT3DAVHOLE Plot Files (IPLOT = 2) 

With the IPLOT = 2 option in Proteus , the XYZ and Q files are written in PLOT3D/WHOLE format. 
With this option, the XYZ file is written using the following Fortran statements: 

write 

write 

$ 

$ 

(nplotx) nl,n2,n3 

Cnplotx) ( <(x(il, 12,13), il=l,nl), 12=1, n2), 13=1, n3), 
((CyCil,i2,i3),il = l,nl),i2=l,n2),i3=l,r»3) , 
(C(z(il,i2,i3),il = l,nl),i2=l ,n2) ,i3-l,n3) 

The Q file is written using the following Fortran statements: 

write 

write 

write 

$ 

(nplot) nl,n2,n3 

Cnplot) machr , aoa , rer , tau( 1 , 1 > 1) 

(nplot) ( ( ( (qplotCil > i2, i3, ivar) , il = l , nl 1 , 12=1 ,n2) , 

13=1 ,n3) ,ivar=l ,5) 

The above WRITE statements for the Q file are executed for each time level written into the file. The Q 
file thus consists of multiple sets of data, each containing the computed results at a single time level. The 
XYZ file is written only once. 16 

The parameters written into the file are defined as follows: 

N1 

Number of grid points N\ in the £ direction. 

N2 

Number of grid points N 2 in the rj direction. 

N3 

Number of grid points A3 in the £ direction. 

X 

Cartesian x coordinate. 

Y 

Cartesian y coordinate. 

Z 

Cartesian z coordinate. 

MACHR 

Reference Mach number, M r — Url^Jy r RTr . This is a real variable. 

AOA 

In PLOT3D, the angle of attack. Set equal to 0. in Proteus. 

RER 

Reference Reynolds number, Re r = p r u r L r l\x r . 

TAU( 1,1,1) 

The time n,i,i. 17 

QPLOT(,„l) 

Static density p. 

QPLOT(,„2) 

x-momentum puM r . 

QPLOT(,„3) 

^-momentum pvM r . 

QPLOT(,„4) 

z-momentum pwM r . 

QPLOT(,„5) 

Total energy per unit volume E r M}. 


16 The current version of PL0T3D does not work for multiple time levels, although future versions might. 

17 Note that with IDTAU = 5 or 6, ? will vary in space, and therefore # ^ 1 , 1 , 1 - 
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All of the parameters written into the XYZ and Q files are nondimensional. However, PLOT3D as- 
sumes that velocity is nondimensionalized by the reference speed of sound a, = (y r RT r Y' 2 , and that energy 
is nondimensionalized by p r a}. In Proteus these variables are nondimensionalized by u r and p r u}. That is 
why M, appears in the definitions of QPLOT(,,,2) through QPLOT(„,5). 

4.2.3 PLOT3D/PLAXES Plot Files (IPLOT = 3) 

With the IPLOT = 3 option in Proteus , the XYZ and Q files are written in PLOT3D/PLANES format. 
With this option, the XYZ file is written using the following Fortran statements: 

write (nplotx) nl,n2,n3 
do 10 i3 = l,n3 

write (nplotx) (CxCil ,i2,i3) ,il = l ,nl) , 12=1 ,n2) , 

$ (Cy(il,i2,i3),il=l,nl),i2=l ,n2) , 

$ (Cz(il,i2,i3),il=l,nl),i2-l,n2) 

10 continue 

The Q file is written using: 

write (nplot) nl,n2,n3 

write (nplot) machr^aoa , rer,tau(l ,1,1) 
do 20 i3 = l,n3 

write (nplot) (((qplotCil, i2, i3,ivar) , il=l ,nl ) , 12=1 , n2) , 

$ ivar=l,5) 

20 continue 

As in the IPLOT = 2 option, the above WRITE statements for the Q file are executed for each time level 
written into the file. The Q file thus consists of multiple sets of data, each containing the computed results 
at a single time level. The XYZ file is written only once. 

4.3 CONVERGENCE HISTORY FILE 

In Section 4.1.6, the convergence history printout is described. The information in this printout is read 
from an unformatted convergence history file that is updated whenever convergence is checked during a 
Proteus calculation. Plotting routines are also available at NASA Lewis to plot any of the convergence 
parameters as a function of time level. 

The file is written in subroutine RESID using the following Fortran statements. At the first time step 
of the run, 


write (nhist) nl ,n2 ,n3,neq,idtau,ictest , ni tavg , ihstag , 

$ iav2e , iav4e , ur , lr , (eps( ieq) , ieq=l , neq) 

Then, at each time level that convergence is being checked, 

write (nhist) i t , tau( 1 , 1 , 1 ) , icheck 

write (nhist) ( chgmax( ieq , 1 ) ,chgavg(ieq) , res!2( ieq , 1 ) , 

$ resavg( ieq, 1 ) , resmax( ieq , 1 ) , lrmaxd , ieq,l) , 

$ lrmax(2 , ieq, 1 ) , lrmax( 3 , ieq , 1 ) ,ieq=l ,neq) 

Finally, again at each time level that convergence is being checked, but only for cases run with explicit ar- 
tificial viscosity, 


write (nhist) ( chgmax( ieq , 1 ) ,chgavg(ieq) , resl2 ( ieq , 2) , 

$ resavgC ieq , 1 ) , resmaxC ieq, 1 ) , lrmax(l ,ieq,2) , 

$ lrmax(2 , ieq, 2) , lrmax(3,ieq , 2) , ieq=l ,neq) 

The parameters written into the file are defined as follows: 

N 1 Number of grid points N\ in the £ direction. 

N2 Number of grid points N 2 in the rj direction. 
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X3 

NEQ 

IDTAU 

ICTEST 

NITAVG 


IHSTAG 

IAV2E 

IAV4E 

UR 

LR 

EPS(IEQ) 

IT 

TAU(1,1,1) 

ICHECK 

CHGMAX(IEQ,1) 

CHGAVG(IEQ) 

RESL2(IEQ,IAVR) 

RESAVG(IEQ,IAVR) 

RESMAX(IEQJAVR) 

LRMAX(IDIR,IEQ,IAVR) 


Number of grid points A 3 in the £ direction. 

The number of coupled governing equations N fq being solved. 
Flag for time step selection method. 

Flag for type of convergence test. 

Number of time steps used in the moving average convergence 
test. 

Flag for constant stagnation enthalpy option. 

Flag for second-order explicit artificial viscosity. 

Flag for fourth-order explicit artificial viscosity. 

Reference velocity Ur . 

Reference length L r . This is a real variable. 

Value e used to determine convergence. 

Current time level n . 

Current value of the time marching parameter t* i i at 

( = *1 = C = 0. 

Results are checked for convergence every ICHECK'th time level. 

Absolute value of the maximum change in the dependent variables 
from time level n — 1 to n. 

Average of the absolute value of the maximum change in the de- 
pendent variables for the last NITAVG time steps. 

The Li norm of the residual at time level n. 

The average absolute value of the residual at time level n. 

The maximum absolute value of the residual at time level n 

The grid indices (ij, k) corresponding to the location of R w „. 


in the above definitions, the subscript IEQ = 1 to N eqf corresponding to the N eq governing equations, 
IAVR = 1 or 2, corresponding to residuals computed without and with the artificial viscosity terms, and 
IDIR = 1, 2, or 3, corresponding to the £, and £ directions. 


4.4 RESTART FILES 


It's often necessary or desirable to run a given case in a series of steps, stopping and restarting between 
each one. This may be done because of limitations in computer resources, or to change an input parameter. 
This capability is provided in Proteus through a restart option. With this option, the computational mesh 
and the computed flow field are written as unformatted output files at the end of one run, saved on a 
magnetic disk or tape, and read in as input files at the beginning of the next run. This process is controlled 
by the input parameters in namelist RSTRT. (See Section 3.1.3). 

The restart files are written and read in subroutine REST. The computational mesh is written using the 
following Fortran statements: 

write Cnrxout) nl,n2,n3 

write (nrxout ) (((x(il,i2,i3),il=l,nl),i2=l,n2),i3=l,n3), 

$ <(Cy(il,i2,i3),il=l,nl),i2=l,n2),i3=l,n3), 

$ (C(z(il ,i2,i3) , il=l ,nl ) , i2=l ,n2) , i3=l ,n3) 

The computed flow field is written using: 

write (nrqout) nl,n2,n3 


60 4.0 Output Description 


Proteus 3-D User's Guide 



write Cnrqout) machr , aoa , rer , tlevel 

write (nrqout) (CCCq C il , 12, 13, ivar) ,11 = 1 >nl) ,i2=l , n2 ) , 

$ i3=l ,n3) ,ivar=l ,5) 

if (iturb .gt. 19) write Cnrqout) (CCCqt ( il , i2 , i3 , ivar ) , 11 = 1 , nl ) , 
$ i2=l ,n2) ,i3=l ,n3) , ivar=l ,2) 

write Cnrqout) ( ( ( CqlC il , i2 , i3, ivar) ,11 = 1 , nl ) , i2=l *n2 ) , 

$ 13=1 ,n3) ,i3=l ,n3) , ivar=l ,5) 

if Citurb .gt. 19) write Cnrqout) C C C CqtlC il , i2 , i3 , ivar ) , il = l ,nl ) , 
$ i2=l ,n2) , i3=l , n3) , ivar=l ,2) 


The parameters written into these files are defined as follows: 


XI 

N2 

N3 

X 

Y 

Z 

MACHR 

AOA 

RER 

TLEVEL 


Number of grid points A\ in the £ direction. 

Number of grid points N 2 in the rj direction. 

Number of grid points Nz in the £ direction. 

Cartesian x coordinate. 

Cartesian y coordinate. 

Cartesian z coordinate. 

Reference Mach number, M, = UrlJy,R T, . This is a real variable. 
Set equal to 0. 

Reference Reynolds number, Re r = p r UrL r lp r . 

The current time level n. 


QCnl) 

QU2) 

QU3) 

QLA) 

Q(n,5) 

QL(,„l-5) 


Static density p at time level n. 
x-momentum puM r at time level n. 
jy-momentum pvM r at time level n. 
z-momentum pwM r at time level n . 

Total energ> T per unit volume EtM} at time level n. 
Same as Q(,„l-5), except at time level n — 1. 


QT(,,,1) Turbulent kinetic energy k at time level n. 

QT(,,,2) Turbulent dissipation rate e at time level n . 

QTL(,,,l-2) Same as QT(„,l-2), except at time level n— 1. 

Note that, except for the QT, QL, and QTL variables, these files have the same format as the XYZ and 
Q files created using the I PLOT = 2 option. These restart files can thus also be used as XYZ and Q files 
for the PLOT3D plotting program. The QT, QL, and QTL variables will not be read by PLOT3D. Note 
also, however, that the reverse is not true. The XYZ and Q files created using the IPLOT = 2 option may 
not be used as restart files, since they do not include the QT, QL, and QTL variables. 18 


18 Actually, if the input parameters IPLT and IPLTA are such that only the final time level is written into the Q file, 
the XYZ and Q files may be used for an "approximate' restart. In this case, Proteus will set QL = Q. 
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5.0 INITIAL CONDITIONS 


Initial conditions are required for the mean flow equations throughout the flow field to start the time 
marching procedure. Although the best choice for initial conditions will be problem-dependent, some 
general comments can be made. First, for unsteady flows they should represent a real flow field. A con- 
verged steady-state solution from a previous run would be a good choice. Second, to minimize the number 
of iterations required for steady flows, the initial conditions should be reasonably close to the expected final 
solution. And third, to minimize problems with starting transients it is important that they represent a 
physically realistic flow. 

If the k-z turbulence model is being used, then initial conditions are also required for k and z throughout 
the flow field. Like the mean flow equations, the best choice for the k-z initial conditions will be highly 
problem-dependent, and will have great influence on the convergence and the stability of the computation. 
It should be noted that k = 0 and z = 0 is merely a trivial solution to the k-z equations (Nichols, 1990), so 
the calculation should not be started with k = z = 0 everywhere in the flow field. The initial profiles for k 
and z should be realistic and smooth, and they should at least obey the boundary conditions (i.e., approach 
zero at solid walls, have zero gradient in the far field, etc.) The initial conditions for the mean flow' prop- 
erties (p f u , v, w, and n t ) should also be realistic for fully turbulent flows. A good place to start a calculation 
with the k-z turbulence model would be a restart run from a converged solution of the same flow field with 
an algebraic turbulence model. 

Initial conditions for the mean flow equations may be supplied by a user- written subroutine, called 
INIT, by a default version of INIT that specifies uniform flow with constant flow properties everywhere in 
the flow field, or by restart mesh and flow field files written during a previous run. 

Initial conditions for the k-z equations may be supplied by a user- written subroutine, called KEINIT, 
by a default version of KEINIT that computes k and z from the mean flow field and the assumption of local 
equilibrium (i.e., production equals dissipation), or by restart mesh and flow field files written during a 
previous k-z calculation. 

5.1 USER-WRITTEN INITIAL CONDITIONS 

5.1.1 Mean Flow Equations 

As stated above, the best choice for initial conditions will be problem-dependent. For this reason 
Proteus does not include a general-purpose routine for setting up initial conditions. Instead, provision is 
made for a user- written subroutine, called INIT, that sets up the initial conditions. 

The time-marching algorithm used in Proteus requires initial conditions for p, u, v, w, £V, p , and TP 
These variables may, of course, be computed from many different combinations of known parameters. To 
make this process reasonably flexible, the user may choose from several combinations exactly which vari- 
ables subroutine INIT will supply. This choice is determined by the input parameter ICVARS in namelist 
FLOW. The following tables list the flow field variables to be supplied by subroutine INIT for the various 
ICVARS options, along with the Fortran variables into which they should be loaded. 20 Remember that the 
initial conditions must be nondimensionalized by the reference conditions listed in Table 3-2. The default 
value for ICVARS is 2. 


19 Fewer variables may actually be needed, depending on the value of the input parameter 1HSTAG. 


20 Note that some input variables, like the Mach number M, are not normally saved as separate Fortran variables. 
To save storage they are to be loaded into existing Fortran variables in INIT. These Fortran variables will later 
be loaded with their normal values. 
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When the energy equation is being solved (IHSTAG = 0), the allowed values are: 


ICVARS Variables Supplied By INIT Fortran Variables 


1 p, pit, pv, pw, Et 

2 p, u, v, w, T 

3 p, u, v, w, 7 

4 p, u, v, w, p 

5 c p , u, v, w, T 

6 p, M, a„ oc„, T 


RHO, U, V, W, ET 
P, U, V, W, T 
RHO, U, V, W, T 
P, U, V, W, RHO 
P, U, V, W, T 
P, U, V, W, T 


When constant stagnation enthalpy is assumed (IHSTAG = 1), the allowed values are: 

ICVARS Variables Supplied By INIT Fortran Variables 

1 p, pu, pv, pw RHO, U, V, W 

2 p, u, v, w P, U, V, W 

3 p, u, v, w RHO, U,V,W 

5 c p , u, v, w P, U, V, W 

6 p, M, a v) P, U, V, W 

In the above tables, c p , a„ and a w represent static pressure coefficient, flow angle in degrees in the x-y 
plane, and flow angle in degrees in the x-z plane, respectively. These parameters are defined as: 

P~Pr 

Cp ~ 2 n 
PPh 


* — 1 V 

a v = tan — 


a w = tan 


w 

u 


The Proteus subroutine INITC will use the variables supplied by subroutine INIT to compute p, u y v, 

and Et, using perfect gas relationships if necessary. From these variables, the pressure and temperature 
will then be recomputed using the equation of state in subroutine EQSTAT, overwriting any values speci- 
fied by the user in INTT. This ensures a consistent set of initial conditions for the time marching procedure. 

Subroutine INTT is called once, immediately after the input has been read, the reference and normalizing 
conditions have been set, and the geometry and mesh have been defined. The user may do anything he or 
she desires in the subroutine, such as reading files, reading additional namelist input, making computations, 
etc. Any of the defined Fortran variables in common blocks may be used. 21 (Of course, they should not 
be changed.) The only requirement is that subroutine INTT return to the calling program (\vhich is INITC) 
the combination of variables specified by ICVARS, defined at every grid point. 

Subroutine INIT is also a convenient place to specify point-by-point boundary condition types and 
values. It's often easier to do this using Fortran coding rather than entering each value into the namelist 
input file. 

5.1.2 k-z Equations 

As with the mean flow equations, the best choice for initial conditions for the k-z equations will be 
problem-dependent. For this reason, Proteus does not have a general-purpose routine for setting up initial 
conditions for k and z. Instead, provision is made for a user-written subroutine, called KEINIT, that sup- 
plies the initial values for k and z . 


21 See Volume 3 for definitions of all the common block variables. 
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Assuming that the variables p, u, v, w, Er , p , and T at time level n are available for the turbulent flow 
field under consideration, the k-z turbulence model in Proteus also requires k y s, and p t values at tune levels 
n and n- 1 before the £-c calculation can be started. These variables may be computed by using many 
different assumptions and techniques, or simply by fitting the available experimental data. 

Subroutine KEINIT is used when a calculation is first started from an initial flow field, or when a restart 
is made from an algebraic turbulence model calculation. It is called once, immediately after the initial 
conditions for the mean flow have been computed. Subroutine KEINIT must supply the initial values for 
the Fortran variables KE, E, \IUT (omit \1UT if the calculation is a restart from a previous algebraic 
turbulence model calculation), KEL, EL, and MUTL at every grid point. Usually, KEL, EL, and MUTL 
can be set equal to KE, E, and MUT, respectively. The default version of KEINIT, described below, is a 
good place to start when developing a user- written version. 

5.2 DEFAULT INITIAL CONDITIONS 


5.2.1 Mean Flow Equations 

A default version of subroutine INIT is built into Proteus that specifies uniform flow with constant flow 
properties everywhere in the flow field. It uses the ICVARS = 2 option (the default value) and reads initial 
flow field values of />, u, v, w, and T from namelist IC. The defaults for these parameter^ are 1.0, 0.0, 0.0, 
0.0, and 1.0, respectively, resulting in an initial flow field with p = p r , u = v = w = 0, and T — T r . If a value 
for ICVARS other than 2 is set in the input, a warning message is generated and ICVARS is reset to 2. 

5.2.2 k-z Equations 

A default version of subroutine KEINIT is built into Proteus that calculates the k and z initial profiles 
using the initial mean flow field values of p, w, v, w, and p, along with the assumption of local equilibrium 
(i.e., production equals dissipation). The turbulent viscosity value is also needed, and is obtained from 
the Baldwin-Lomax algebraic model. The assumption of local equilibrium has been found to give good 
initial values for the k profile near a solid wall. In particular, the location and magnitude of the first peak 
in the k profile is well predicted. Of course, k and z will depend heavily on the initial velocity and turbulent 
viscosity, and unrealistic initial profiles for k and z might result from bad initial velocity and turbulent 
viscosity profiles. The k and z values calculated using this subroutine can be extremely small in regions of 
nearly inviscid, freestream flow. Therefore, it is up to the user to make sure that the freestream values of 
k and z are appropriate for the problem under consideration, perhaps by modifying the value of the input 
parameter CKMIN. 

5.3 RESTART INITIAL CONDITIONS 

If restart mesh and flow field files were created during a previous run by setting IREST > 0 in namelist 
RSTRT, these files can be used to continue the calculation from the point where the previous run stopped. 
In this case no subroutine INIT is needed. The restart case is run by linking the existing restart mesh and 
flow field files to Fortran units NRXIN and NRQIN, respectively, and setting IREST — 2 or 3 in the input. 
New restart files will also be written to units NRXOUT and NRQOUT. 

As mentioned above, subroutine KEINIT is used when the k-z computation is done for the first time 
(IREST < 1), or when restarting from a previous computation that used an algebraic turbulence model 
(IREST - 3). However, if the previous calculation used the k-z turbulence model (IREST = 2), then the 
restart files also contain the information needed to continue the k-z calculation. In this case, subroutine 
KEINIT is not needed. 

When a restart case is being run, the usual namelist input described in Section 3.1 must still be read in. 
The following input parameters must have the same values as during the original run: IUNITS, IHSTAG, 
ILAMV, LR, UR, MACHR, TR, RHOR, MUR, RER, KTR, PRLR, GAMR, RG, HSTAGR, Nl, N2, 
N3, IPACK, and SQ. The remaining input parameters either may be changed during a restart, or do not 
apply to a restart case. Note, however, that for many of the input parameters, such as those^ specifying the 
boundary conditions, changing values during a restart may result in temporary re-starting transients or 
even cause the calculation to blow 7 up. 
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6.0 RESOURCE REQUIREMENTS 


Proteus was developed on the Cray X-MP and Y-MP computers at NASA Lewis Research Center. 
Changes that may be necessary when porting the code to another system are described in Section 6.1. 
Sections 6.2 and 6.3 discuss the memory and CPU time required to run the code. The values presented in 
these sections are for version 1.0 of the 3-D Proteus code, and were derived from tests run on the NASA 
Lewis Cray Y-MP in March, 1992. At that time the Cray was running UNICOS Version 6.0 and CFT77 
5.0.2. UNICOS and CFT77 are described in the UNICOS User Commands Reference Manual (Cray Re- 
search, Inc., 1991), and in CF77 Compiling System, Volume 1: Fortran Reference Manual (Cray Research, 
Inc., 1990), respectively. Section 6.4 describes the size and format of the various input and output files used 
in the code. 

6.1 COMPUTER 


Proteus should be transportable to other computer systems with minimal changes. With just three 
known exceptions, the code is written entirely in ANSI standard Fortran 77. The first exception is the use 
of namelist input. With namelist input, it's relatively easy to create and/or modify input files, to read the 
resulting files, and to program default values. Since most Fortran compilers allow namelist input, its use 
is not considered a serious problem. The second exception is the use of + CALL statements to include 
COMDECK's, which contain the labeled common blocks, in most of the subprograms. This is a Cray 
UPDATE feature, and therefore the source code must be processed by UPDATE to create a file that can 
be compiled. 22 UPDATE is described in the UPDATE Reference Manual (Cray Research, Inc., 1988). Since 
using the +CALL statements results in cleaner, more readable code, and since many computer systems have 
an analogous feature, the + CALL statements were left in the program. The third exception is the use of 
lowercase alphabetic characters in the Fortran source code. This makes the code easier to read, and is a 
common extension to Fortran 77. 

Several library subroutines are called by Proteus. SGEFA and SGESL are Cray versions of LINPACK 
routines. SASUM and SNRM2 are Cray Basic Linear Algebra Subprograms (BLAS). GATHER is a Cray 
Linear Algebra routine. ISAMAX, ISAMIN, ISRCHEQ, ISRCHFGT, ISR^HFLT, and WHENFLT are 
Cray search routines. TREMAIN is a Cray Fortran library routine. All of these routines are described in 
detail in Section 4.0 of Volume 3. In addition, SGEFA and SGESL are described in Volume 3: UNICOS 
Math and Scientific Library Reference Manual (Cray Research, Inc., 1989b) and by Dongarra, Moler, 
Bunch, and Stewart (1979); SASUM, SNRM2, GATHER, ISAMAX, ISAMIN, ISRCHEQ, ISRCHFGT, 
ISRCHFLT, and WHENFLT are described in Volume 3: UNICOS Math and Scientific Library Reference 
Manual (Cray Research, Inc., 1989b); and TREMAIN is described in Volume I: UNICOS Fortran Library 
Reference Manual (Cray Research, Inc., 1989a). These or similar routines may be available on other sys- 
tems. If not, equivalent routines will have to be coded. 

6.2 MEMORY 


The sizes of the dimensioned arrays in Proteus, and hence the amount of memory required to run the 
program, are set using Fortran parameter statements. These parameters are set in COMDECK PARAMS1. 
Larger or smaller dimensions may be set for the entire program simply by changing the appropriate pa- 
rameter, and then recompiling the entire program. The parameters are defined as follows: 

NIP Maximum number of grid points in the c, direction. The current value is 21. 

N2P Maximum number of grid points in the r\ direction. The current value is 21. 

N3P Maximum number of grid points in the £ direction. The current value is 21. 


22 See the example in Section 8.1. PRECEDING PAGE BLANK NOT FfLNtED 
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NMAXP 

N'TOTP 

NDIAGP 

NPNTP 

NEQP 


NEQPM 

NAMAX 

NBC 

NTP 

NTSEQP 


Maximum of NIP, N2P, and N3P. 

Total storage required for a single three-dimensional array (i.e., 

NTP x N2P x N3P). 

Number of diagonals containing interior points in a plane (i.e., 

NTP + N2P-5). 

Number of interior points in a %-rj plane (i.e., (NTP — 2)(N2P — 2)). 

Number of coupled equations allowed. The number of equations to be solved 
depends on the value of the input parameter IHSTAG, as shown in Table 3-4 in 
Section 3.0. In the Proteus code, NEQP is initially set equal to 5. Note that 
Proteus will still work if NEQP is larger than the value of NEQ shown in Table 
3-4. 

Maximum number of coupled equations available. This is the largest permissible 
value for NEQP- The current value is 5. 

Maximum number of time steps allowed in the moving average convergence test 
(the ICTEST = 2 option). The current value is 10. 

Number of boundary conditions per equation. The current value is 2. 

Maximum number of entries in the table of time-dependent boundary condition 
values. The current value is 10. 

Maximum number of time step sequences in the time step sequencing option. The 
current value is 10. 


The total amount of memory in computer words required to run Proteus 3-D, compiled using CFT77, 
is listed in the following table for various combinations of the dimensioning parameters NTP, N2P, N3P, 
and NEQP. 23 On the Cray X-MP and Y-MP, each word is 64 bits long. 


NIP 

N2P 

N3P 

MEMORY 

NEQP =4 

NEQP = 5 

21 

21 

21 

901,632 

913,920 

51 

51 

51 

7,016,960 

7,089,664 

101 

51 

51 

13,721,088 

14,00,, 784 


6.3 CPU TIME 


Compilation of Proteus 3-D using the CFT77 compiler requires about 235 seconds of CPU time. The 
CPU time required for execution will depend on several factors, such as the turbulence model being used, 
and whether or not the energy equation is being solved. For the test case in Section 9.0, the CPU time 
ranged from about 4.8 x 10“ 5 to 5.7 x 10“ 5 seconds per grid point per time step. 

6.4 INPUT/OUTPUT FILES 

Several files are used in Proteus for various types of input and output. The contents of these files have 
been described previously in Sections 3.0 and 4.0. This section describes the characteristics of the files 
themselves. The fdes are identified by the Fortran variable representing the unit number. The unit numbers 
have default values, but all of them except NTN may be read in by the user. 

Table 6- 1 lists the files used in Proteus, giving the default unit number, briefly describing the contents 
of the file, and indicating when it is used. Table 6-2 summarizes the computational resources needed for 
each file. In this table, the record length is specified in bytes for units NTN and NOUT, and in computer 
words for the remaining units. The total file size is specified in printed pages for unit NOUT, and in 


23 It should be noted that version 5.0.2 of the CFT77 compiler itself requires just under 2,400,000 words. 
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computer words for the remaining units. Several symbols and Fortran variables are used in Table 6-2, and 
are defined as: 


N m 

A rurb 

K, 


Art. A a 

N(, A r „ N ( 


M. Ay, A 3 

ICHECK 

IPLOT 

NG1, NG2, NG3 


0 if explicit artificial viscosity is not being used, 1 if it is. 

0 if the input parameter ITURB = 0 or 1, 1 if ITURB = 20. 

The number of coupled equations being solved. 

The number of time levels written into the plot file. This is determined 
by the input parameters IPLT and IPLTA. 

The time level at the beginning and at the end of the calculation. 

The number of grid points in the £, »j, and C directions at which output is 
being printed. This is determined by the input parameters IPRT1, IPRT2, 
IPRT3, IPRT1A, IPRT2A, and IPRT3A. 

The number of grid points in the <*, »j, and C directions. 

Input parameter specifying frequency for checking convergence. 

Input flag specifying type of plot file being written. 

Number of points in the £ , V. and £ directions in the coordinate system file. 


The typical record length and total size values listed in the table were computed assuming N„= 1, 
Nturb = o, N eq = 5, N, = 1, Art = 1, Na = 2000, N { = A, = A' { = 21, A, = N 2 = A 3 = 51, ICHECK = 10, and 
NG1 = NG2 = NG3 = 51. 


Proteus 3-D User's Guide 


6.0 Resource Requirements 69 


TABLE 6-1. - I/O FILE TYPES 


UNIT 

DEFAULT 
UNIT NO. 

RECORD 

FORMAT 

CONTENTS 

WHEN USED 

NIN 

5 

Formatted 

Standard input 

Always 

NOUT 

6 

Formatted 

Standard output 

Always 

NGRID 

7 

Unformatted 

Coordinate system input 

NGEOM = 10 

NPLOTX 

8 

Unformatted 

PLOT3D XYZ file output 

I PLOT = 2, 3 

NPLOT 

9 

Unformatted 

CONTOUR plot file or 
PUOT3D Q file output 

IPLOT * 0 

NHIST 

10 

Unformatted 

Convergence history output 

Always 

NRQIN 

11 

Unformatted 

Restart flow field input 

IREST = 2, 3 

NRQOUT 

12 

Unformatted 

Restart flow field output 

IREST = I, 2, 3 

NRXIN 

13 

Unformatted 

Restart computational 
coordinates input 

IREST = 2 , 3 

NRXOUT 

14 

Unformatted 

Restart computational 
coordinates output 

IREST = 1, 2, 3 


70 6.0 Resource Requirements 


Proteus 3-D User's Guide 










TABLE 6-2. - I/O FILE SIZES 


UNIT 

MAXIMUM 
RECORD LENGTH 4 

TYPICAL 

MAXIMUM 

RECORD 

LENGTH 4 

TOTAL SIZE* 5 

TYPICAL 

TOTAL 

SIZE b 

NIN 

80 


Variable 

Variable 

NOUT 

132 


132 

=; (2.V„ + 3) pages + 

13 pages + 





ty{(Y, + 3)[(A<- 1)/10+ 1] 

42 pages per 





-f- 2}/55 pages per variable 

variable per 





per time level 

time level 

NGRID 

3(NG1)(NG2)(NG3) 

397,953 

3(NG1)(NG2)(NG3) + 2 

397,955 

NPLOTX 

IPLOT 


I 

IPLOT 




2 

3N l N 2 N l 

■ 

2 

3NiN 2 N 3 + 3 

397,956 


3 

3 A’,JV 2 


3 

3AvV 2 N 3 + 3 

397,956 

NPLOT 

IPLOT 



IPLOT 




1 

14;V 3 

714 


(14YiY 2 A 3 + 23) N, 

1,857,137 


2 

SNuViNi 

663,255 


(5A)A 2 A 3 + 1)N, 

663,262 


3 

5NiN 2 

13,005 


(5NiN 2 Ni + 7 )N, 

663,262 

NHIST 

8;V„ 


40 

[(A„ + l)(8Y e ,) + 3] ■ 

8,617 





(N,2 - N n + 1)/ICHECK 






+ N" + 12 


NRQIN, 

5NiN 2 N 3 


663,255 

lOY.AW, + 7 

1,326,517 

NRQOUT 







NRXIN, 

3NiN 2 Nj 


397,953 

3A,Y 2 A 3 + 3 

397,956 

NRXOUT 








a In bytes for units NIN and NOUT 7 and in computer words for the remaining units. 
b In pages for units NIN and NOUT t and in computer words for the remaining units. 
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7.0 DIAGNOSTIC MESSAGES 


Various diagnostic messages may be printed by Proteus as part of its standard output. -4 Most of these 
concern inconsistent or invalid input, although some describe problems encountered during the calculation 
itself. Two types of messages may appear - errors and warnings. Error messages are preceded by the 
characters **** ERROR, and indicate either a serious problem or an input error that Proteus cannot resolve. 
Warning messages are preceded by the characters **** WARNING, and indicate either a potential problem 
or a non-standard combination of input parameters. Errors cause the calculation to stop, while warnings 
do not. 

The various error and warning messages are listed and explained in the following two subsections. Italic 
letters, like value, are used to indicate variable values. The subprogram in which the message is printed is 
given in parentheses at the end of the each explanation. 

7.1 ERROR MESSAGES 


Both machr and ur specified, 
mac hr = value, ur = value 

Either the reference Mach number or velocity may be specified in namelist FLOW, but not both. 
(INPUT) 

Both prlr and ktr specified, 
prlr = value, ktr = value 

Either the reference laminar Prandtl number or thermal conductivity may be specified in namelist 
FLOW, but not both. (INPUT) 

Both rer and mur specified. 

rer = value, mur = value 

Either the reference Reynolds number or viscosity may be specified in namelist FLOW, but not 
both. (INPUT) 

Coordinate system file has ngl, ng2, and/or ng3 > max allowed, 
ngl = value, ng2 = value, ng3 = value 

A coordinate system file has been read in, using the NGEOM = 10 option, with more grid points 
than allowed. The maximum allowed values of NG1, NG2, and NG3 are the values of the di- 
mensioning parameters NIP, N2P, and N3P, respectively. (GEOM) 

Grid transf ormation Jacobian changes sign or = 0 . 

The nonorthogonal grid transformation Jacobian J must be either everywhere positive or every- 
where negative. This error indicates that the computational mesh contains crossed or coincident 
grid lines. The error message is followed by a printout of the Cartesian coordinates, the Jacobian, 
and the metric coefficients. (METS) 

Illegal option for computational coordinates requested, 
ngeom = value 

An illegal value of NGEOM has been specified in namelist GMTRY. The legal values are 1,2, 
and 10, and are described in Section 3.1.5. (GEOM) 


24 The diagnostic messages described in this section are generated by the Proieus code itself, and appear as part of 
the standard output. Any computer system error messages due to floating-point errors, etc., are, of course, 
system-dependent. On UNIX-based systems, system errors will normally appear in the standard error file. 
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Illegal plot file option requested. 

ip lot = value 

An illegal value of IPLOT has been specified in namelist 10. The legal values are 0, 1, 2, and 3, 
and are described in Section 3.1.4. (PLOT) 

Illegal thin-layer option, 
ithin lidir) = value 

An illegal value of YTHlS(idir) has been specified in namelist FLOW. Here idir = 1, 2, or 3, 
corresponding to the £, 77 , and £ directions, respectively. The legal values are 0 and 1 , and are de- 
scribed in Section 3.1.6. (INPUT) 

Illegal time step selection option requested, 
idtau = value 

An illegal value of IDTAU has been specified in namelist TIME. The legal values are 1 to 9, and 
are described in Section 3.1.9. (TIMSTP) 

Illegal value for icvars. 
icvars = value 

An illegal value of ICVARS has been specified in namelist FLOW. The legal values are 1 to 6 , and 
are described in Section 3.1.6. (INITC) 

Illegal value for ieuler. 
ieuler = value 

An illegal value of IEULER has been specified in namelist FLOW. The legal values axe 0 and 1, 
and are described in Section 3.1.6. (INPUT) 

Illegal value for ihstag. 
ihstag = value 

An illegal value of IHSTAG has been specified in namelist FLOW. The legal values are 0 and 1, 
and are described in Section 3.1.6. (INPUT) 

Illegal value for ilamv. 
ilamv = value 

An illegal value of ILAMV has been specified in namelist FLOW. The *egal values are 0 and 1, 
and are described in Section 3.1.6. (FTEMP) 

Illegal value for ildamp. 
ildamp = value 

An illegal value of ILDAMP has been specified in namelist TURB. The legal values are 0 and 1, 
and are described in Section 3.1.10. (INPUT) 

Illegal value for inner, 
inner = value 

An illegal value of INNER has been specified in namelist TURB. The legal values are 1 and 2 , and 
are described in Section 3.1.10. (INPUT) 

Illegal value for irest. 
irest = value 

An illegal value of IREST has been specified in namelist RSTRT. The legal values are 0, I, 2, and 
3, and are described in Section 3.1.3. (INPUT) 

Illegal value for iunits. 
i units = value 

An illegal value of IUNITS has been specified in namelist 10. The legal values are 0 and 1 , and 
are described in Section 3.1.4. (INPUT) 
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Illegal value for lwalll. 

IwallK j ,k, ibound) = value 

Illegal value for lwall2. 

Iwall2 (i,k, ibound) = value 

Illegal value for lwall3. 

Iwall3(/ ,j Abound) = value 

An illegal value of LWALL1, LWALL2, or LWALL3 has been specified in namelist TURB. Here 
ij y and k are the indices in the £, rj, and £ directions, and ibound = 1 or 2, corresponding to the 
£ = 0 or 1 surface, the ^ = 0 or 1 surface, or the £ = 0 or 1 surface. The legal values are 0 and 1, 
and are described in Section 3.1.10. (INPUT) 

Improper specification of thin*layer option, 
ithin = value value value 

The thin-layer option may be used in one direction only. Thus, one element of the ITHIN array 
must be set equal to 1, and the remaining two must be 0. (INPUT) 

Invalid boundary condition type requested. 
jbclUeq Abound) or ibcl (j ,k , ieq Abound) = value 

Invalid boundary condition type requested. 
jbcZUeq Abound) or ibc2( i ,k ,ieq , ibound) = value 

Invalid boundary condition type requested. 
jbc3( ieq Abound) or ibc3 ii ,j ,ieq , ibound) = value 

These messages result from an invalid boundary condition type being specified in namelist BC for 
the £, rj f and/or £ direction. Here ieq is the boundary condition equation number; ibound = 1 or 
2, corresponding to the £ = 0 or 1 surface, the rj = 0 or I surface, or the £ = 0 or 1 surface; and i, 
/, and k are the indices in the £, yj, and £ directions. The valid boundary conditions are listed in 
Table 3-6. (BCDENS, BCF, BCGEN, BCNVEL, BCPRES, BCQ, BCTEMP, BCUVEL, 
BCWEL, BCWVEL, BC1VEL, BC2VEL, BC3VEL) 

Invalid boundary type requested, 
kbcl abound) = value 

Invalid boundary type requested. 
kbc2( ibound) = value 

Invalid boundary type requested. 
kbc3 abound) - value 

These messages result from an invalid boundary type being specified in namelist BC, for the £, y \ , 
and/or £ direction, when the KBC meta flags are used. Here ibound — 1 or 2, corresponding to 
the £ = 0 or I surface, the rj = 0 or 1 surface, or the £ = 0 or 1 surface. The valid boundary types 
are listed in Section 3.1.7. (BCSET) 

Invalid grid packing location for Roberts formula. 

sq C idir , 1 ) = value 

An invalid grid packing location, given by the value of SQ(IDIR,1) in namelist NUM, has been 
specified. Here idir = 1, 2, or 3, corresponding to the £, yj, and £ directions, respectively. The valid 
values are 0.0, 0.5, and 1.0, and are described in Section 3.1.8. (INPUT) 

Invalid grid packing option. 

ip a ck (idir) = value 

An invalid grid packing option, given by the value of IPACK(IDIR) in namelist NUM, has been 
specified. Here idir = 1, 2, or 3, corresponding to the £, rj, and £ directions, respectively. The valid 
values are 0 and 1, and are described in Section 3.1.8. (PAK) 
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Invalid grid packing parameter for Roberts formula. 
sq(idir,2) = value 

An invalid grid packing parameter, given by the value of SQ(idir,2) in namelist NUM, has been 
specified. Here idir = 1, 2, or 3, corresponding to the r\, and C directions, respectively. The valid 
values are > 1, and are described in Section 3.1.8. (INPUT) 

Invalid time step selection method for time step sequencing option, 
idtau = value , ntseq = value 

A time step selection option that adjusts At as the solution proceeds has been specified in namelist 
TIME in conjunction with the time step sequencing option. If the time step sequencing option is 
being used, IDTAU must be 1, 3, or 5. (INPUT) 

Invalid type of artificial viscosity requested. 
iav4e, iav2e, iav2i = value value value 

An invalid type of artificial viscosity has been specified in namelist NUM. The valid values are 0, 
1, and 2 for IAV4E and IAV2E, and 0 and 1 for IAV2I, and are described in Section 3.1.8. (IN- 
PUT) 

Invalid type of unsteadiness for boundary condition requested. 

j t b c 1 C ieq , ibound ) = value 

Invalid type of unsteadiness for boundary condition requested. 

5±bc2Ueq , ibound) = value 

Invalid type of unsteadiness for boundary condition requested. 

5 tbcZUeq , ibound) = value 

These messages result from an invalid type of unsteadiness being specified in namelist BC for the 
boundary conditions in the £, and/or ( direction. Here ieq is the boundary condition equation 
number, and ibound = 1 or 2, corresponding to the £ = 0 or 1 surface, the rj — 0 or 1 surface, or 
the £ = 0 or 1 surface. The valid values for JTBC1, JTBC2, and JTBC3 are 0, 1, and 2, and are 
described in Section 3.1.7. (TBC) 

Mesh size requested > max allowed. 

nl = value , n2 = value, n3 = value 

More grid points have been requested in namelist NUM than are allowed. For non-periodic 
boundary conditions, the maximum allowed values of Nl, N2, and N3 are the values of the di- 
mensioning parameters NIP, N2P, and N3P, respectively. For spatially periodic boundary condi- 
tions, the maximum values are NIP — 1, N2P — 1, and N3P — 1. (INPUT) 

More time step sequences requested than allowed, 
ntseq = value 

For the time step sequencing option, the number of time step sequences, specified in namelist 
TIME, cannot exceed the value of the dimensioning parameter NTSEQP. (INPUT) 

neqp not large enough, 
neqp, ihstag = value value 

The value of the dimensioning parameter NEQP, which sets the maximum number of coupled 
equations that can be solved, is not large enough. The number of equations to be solved depends 
on the value of the input parameter IHSTAG. See Section 6.2 and Table 3-4 in Section 3.0. 
(INPUT) 

Non-existent turbulence model requested, 
iturb = value 

A non-zero value for ITURB has been specified in namelist FLOW that does not correspond to 
one of the turbulence models currently available in Proteus . The only valid non-zero values for 
ITURB are I and 20, corresponding to the algebraic Baldwin-Lomax and the Chien two -equation 
turbulence models, respectively. (INITC, MAIN) 
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Non-positive pressure and/or temperature at time level n 

II 12 13 P T 

value value value value value 

During the solution, a non-positive value for pressure and/or temperature has been computed in 
subroutine EQSTAT. Up to 50 values will be printed. These values, of course, are non-physical 
and indicate a failure of the solution. Although the calculation will stop, the standard output and 
plot file will include this time level, if that is consistent with the IPRT and IPLT type parame- 
ters in namelist 10. The restart files will not be written. This failure may be caused by bad initial 
or boundary conditions, or by too large a time step. (INITC, MAIN) 

Number of values in unsteady boundary condition table > max allowed, 
ntbc = value 

For unsteady boundary conditions, the number of values in the tables of GTBC1, GTBC2, and/or 
GTBC3 vs. NTBCA, specified in namelist TIME, cannot exceed the value of the dimensioning 
parameter NTP. (INPUT) 

Singular block matrix for b. c. at lower boundary* sweep n 

Singular block matrix for b. c. at upper boundary, sweep n 

When boundary conditions are specified using the JBC and/or IBC input parameters, zero values 
may appear on the diagonal of the block tridiagonal coefficient matrix. Subroutine FILTER at- 
tempts to rearrange the elements of the boundary condition block submatrices to eliminate any of 
these zero values. These messages indicate it was unable to do so for the boundary and sweep in- 
dicated. This means the diagonal submatrix B is singular, which in turn means the specified 
boundary conditions are not independent of one another. (FILTER) 

7.2 WARNING MESSAGES 

chgmax > 0.15, cfl cut in half. 

chgmax > 0.15, dtau cut in half. 

With the IDTAU = 2, 4, and 6 options, the time step is adjusted gradually as the solution proceeds 
based on the absolute value of the maximum change in the dependent variables. One of these 
messages may occur if the solution changes very rapidly. (The first message applies to the 
IDTAU = 2 and 6 options, and the second to the IDTAU = 4 option.) Under these conditions the 
time step is arbitrarily cut in half, and the solution continues. This may stabilize the calculation, 
but consideration should be given to rerunning the problem with a smaller time step, especially for 
unsteady flows. (TIMSTP) 
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Conflicting boundary conditions specified. 

kbcliibound ) = value, jbcl (ieq, ibound) = value, ibcl( j ,k,ieq, (bound) = value 

kbcliibound ) = value will be used. 

Conflicting boundary conditions specified. 

kbc2( ibound) = value, jbc2 (ieq , ibound) - value, ib c2(i,k, ieq , ibound) = value 

kbc2 (ibound) = value will be used. 

Conflicting boundary conditions specified. 

kbc3 Ubound) = value, jbcZ( ieq , ibound) = value, ibc3( i ,j , ieq , ibound) = value 
kbcZdbound ) = va/ue will be used. 

These messages indicate that a bound ary type was specified using the KB C parameter, and that 
individual boundary conditions were also specified for that boundary using the JBC and/or IBC 
parameters, for one of the £, rj, and/or £ boundaries. For a given boundary, the boundary type 
may be specified using the KBC parameter, or individual boundary conditions may be specified 
using the JBC and/or IBC parameters. Here ibound = 1 or 2, corresponding to the £ = 0 or 1 
surface, the rf = 0 or 1 surface, or the £ = 0 or 1 surface; ieq is the boundary condition equation 

number; and /, y, and k are the indices in the £, r\ ) and £ directions. Boundary' conditions will be 

set using the KBC parameter, and the calculation will continue. See the discussion of boundary 
condition input in Section 3.1.7. (INPUT) 

Conflicting boundary conditions specified, 
jbcl (ieq, ibound) * value, ibclC j ,k , ieq , ibound) = value 
jbcl (ieq, ibound) = value will be used. 

Conflicting boundary conditions specified. 
jbc2 (ieq , ibound) = value, ibc2( i ,k , ieq Abounds - value 
jbc2( ieq, ibound) = value will be used. 

Conflicting boundary conditions specified. 
jbc3 (ieq , ibound) - value, ibc3(/ ,j , ieq , ibound) = value 
jbc3 (ieq , ibound) = value will be used. 

These messages indicate that both surface and point-by-point boundary conditions were specified 
for the same boundary condition equation number for one of the £, and/or £ boundaries. Each 
boundary condition on each boundary may be specified for the entire surface using the JBC and 
GBC parameters, or on a point-by-point basis using the IBC and FBC A irameters, but not both. 
Here ieq is the boundary condition equation number; ibound — 1 or 2, corresponding to the 
£ s= 0 or 1 surface, the y\ = 0 or 1 surface, or the £ *= 0 or 1 surface; and /, j, and k are the indices 
in the £, */, and £ directions. A likely cause of this error message is specifying point-by-point 
boundary conditions without setting the appropriate JBC parameter equal to — 1. The boundary 
condition will be set using the JBC parameter, and the calculation will continue. See the discussion 
of boundary' condition input in Section 3.1.7. (INPUT) 

icvars reset to 2 for default init. 

The default version of subroutine INIT is set up assuming ICVARS = 2. If another value of 
ICVARS is read in, it is automatically reset to 2 and the calculation will continue. (INIT) 

Illegal convergence testing option requested, 
ictest = value, reset to ictest = 3 

An illegal value of ICTEST has been specified in namelist TIME. ICTEST will be reset to 3, 
corresponding to convergence based on the L 2 norm of the residual, and the calculation will con- 
tinue. The legal values are 1 to 5, and are described in Section 3.1.9. (CONV) 

Illegal output requested. 
ivoutC n) = value 

An illegal value of IVOUT has been specified in namelist 10. It will be ignored and the calculation 
will continue. The legal values of IVOUT are listed in Table 3-3. (OUTPUT) 
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Implicit & nonlinear explicit artificial viscosity both requested. 
iav2i, iav2e, iav4e = value value value 
cavs2i - value value value value value 

Normally, Jameson's nonlinear artificial viscosity, specified by setting IAV2E and IAV4E = 2, is 
explicit only. This message is printed when implicit artificial viscosity is requested at the same time. 
Proteus will assume that you know what you are doing and the implicit artificial viscosity will be 
included. (INPUT) 

Invalid debug option specified. 
idebug(Z) 

An invalid debug option, number /, has been specified in namelist 10. The valid IDEBUG options 
are 1 to 7, and are described in Section 3.1.4. (INPUT) 

Non-standard time difference centering requested. 

the = value value 

thx = value value value 

thy = value value value 

thz = value value value 

the = value value value 

The Beam-Warming type time differencing used in Proteus includes three standard implicit schemes 
- Euler, trapezoidal, and three -point backward. This message indicates that a combination of time 
differencing parameters 6 1 , 0 2 , and # 3 has been specified for at least one of the governing equations 
that does not correspond to any of the three standard schemes. Proteus will assume that you know 
what you are doing and the specified time differencing parameters will be used. (INPUT) 

Spatially varying time step requested with time-accurate differencing scheme. 

idtau = value 

the = value value 

thx = value value value 

thy = value value value 

thz = value value value 

the = value value value 

For steady flows, a spatially varying time step may be used to enhan. *. convergence, and first-order 
time differencing is recommended. Using a second-order time-accurate differencing scheme should 
not give wrong results, but is wasteful. For unsteady flows, second-order time-accurate differencing 
should be used, but only a global (i.e., constant in space) time step makes sense. (INPUT) 

Time level may fall outside range of input table for unsteady b. c. 
itbeg = value , itend = value, ntbca(l) = value, ntbcaC/tf&c) = value 

General unsteady boundary' conditions are being used, and the time levels in the input table of 
GTBC1, GTBC2, and/or GTBC3 vs. NTBCA do not cover the time levels that will occur during 
the time marching loop. Here itbeg and itend are the first and last time levels in the time 
marching loop, and ntbc is the value of the input parameter NTBC. If the time level 
/z<NTBCA(l), the boundary' condition value will be set equal to the first value in the table. 
Similarly, if n > NTBCA(NTBC), the value will be set equal to the last value in the table. (TBC) 

7.3 OTHER MESSAGES 


Converged solution at time level n 

The solution has converged, according to the criteria specified by the parameters ICTEST and EPS, 
at the indicated time level. The calculation will stop, and the standard output and plot file will 
include this time level, if that is consistent with the "IPRT" and "IPLT* parameters in namelist 10. 
The restart files will be written. (MAIN) 
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Not yet converged. 

The solution has not converged, according to the criteria specified by the parameters ICTEST and 
EPS, within the number of"time steps specified by NTIME. The calculation will stop, and the 
standard output and plot file will include this time level, if that is consistent with the TPRT" and 
TPLT" parameters in namelist 10. The restart files will be written. (MAIN) 

Sorry, ran out of cpu time at time level n 

The calculation was stopped at the indicated time level because the job w’as almost out of CPU 
time. The standard output and plot file will include this time level, if that is consistent with the 
TPRT" and TPLT 7 parameters in namelist IO. The restart files will be written. (MAIN) 
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8.0 JOB CONTROL LANGUAGE 


At NASA Lewis, Proteus is currently being run on the Cray X-MP and Y-MP computers, with 
UNICOS 6.0 as the operating system. In this section several general examples are presented showing the 
UNICOS job control language (JCL) that may be used as starting points when setting up specific cases. 25 
The individual UNICOS commands are described in detail in the UMCO.S User Commands Reference 
Manual, (Cray Research, Inc., 1991). These examples are written for the Bourne shell. Some changes may 
be necessary if the C shell is being used. It is assumed that the user is familiar with the procedures used at 
their site to submit jobs to the Cray, and to receive and process the output files. 

Each example is given with reference line numbers, which are not part of the actual JCL statement, 
followed by a line-by-line explanation. Note that in UNICOS, the case (upper or lower) of the letters in 
commands and arguments is significant. In this respect, the examples should be followed exactly. 

8.1 COMPILING THE CODE 


In this example, the Proteus code is compiled on the Cray. It is assumed that the source code is in the 
proper format for the Cray UPDATE facility, as described in the UPDATE Reference Manual (Cray Re- 
search, Inc., 1988). 


1 . 

# QSUB - 1M 2 . OmW -IT 

300 

2. 

# QSUB -mb -me 


3. 

# QSUB -o temp.out 


4. 

# QSUB -eo 


5. 

# QSUB -r EXAMPLE 


6. 

# QSUB -s /bin/sh 


7. 

set -x 


8. 

update ~i p3dl0.s -n 

p3dl0.u -c p3dl0 -f 

9. 

cft77 -b p3dl0 . o -d 

pq -o aggress p3d!0.f 

10. 

rm p3dl0.f 



Lines 1 through 6 are actually Network Queueing System (NQS) commands, not UNICOS commands. 
They must appear first in the runstream, and begin with a # sign. This first line sets the memory and CPU 
time limits for the job at 2.0 million words and 300 seconds. 

Line 2 tells the Cray to send email when the job starts and finishes. 

Line 3 tells the Cray to store the standard output in the file temp.out. 

Line 4 causes any system error messages to be included in the standard output file. 

Line 5 gives the name of the job as EXAMPLE. 

Line 6 causes the Boume shell to be used for this job. 

Line 7 causes your UNICOS commands to be printed as part of the output runstream. 

Line 8 uses the Cray UPDATE facility to create the files, p3di0.f which contains the complete compilable 
Fortran code for Proteus , and p3d!0.u, , which contains the Proteus update program library. The update 
command uses as input the file p3dJ0.s . 

Line 9 compiles the program p3d!0.f, storing the object code in the file p3dl0.o. 

Line 10 removes the named file. 


75 See Section 9.0 for specific examples of actual cases. 
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8.2 RUNNING THE MASTER FILE 


The simplest way to run Proteus is shown in this example. The existing master file is being used, 
without making any changes. 

1. # QSUB -1M 4 . OmW -IT 3600 

2. # QSUB -mb -me 

3. # QSUB -o temp. out 

4. # QSUB -eo 

5. # QSUB -r EXAMPLE 

6. # QSUB -s /bin/sh 

7. set -x 

8 . d a 

9. touch plotx 

10. touch plotq 

11. touch chist 

12. touch rq2 

13. touch rx2 

14. In plotx fort. 8 

15. In plotq fort. 9 

16. In chist fort. 10 

17. In rq2 fort. 12 

18. In rx2 fort. 14 

19. In coords fort. 7 

20. segldr -o p3dl0.ex p3dl0.o 

21 . p3dl0 .ex << EOD 

Standard Proteus input file goes here 

22. EOD 

23. rm p3dl0.ex 

24. ja -cslt 

Lines 1 through 6 are actually Network Queueing System (NQS) commands, not UNICOS commands. 
They must appear first in the file, and begin with a # sign. This first line sets the memory and CPU time 
limits for the job at 4.0 milli on words and 3600 seconds. 

Line 2 tells the Cray to send email when the job starts and finishes. 

Line 3 tells the Cray to store the standard output in the file temp.out. 

ijne 4 causes any system error messages to be included in the standard output fik. 

Line 5 gives the name of the job as EXAMPLE. 

Line 6 causes the Bourne shell to be used for this job. 

Line 7 causes your UNICOS commands to be printed as part of the output runstream. 

Line 8 tells the Cray to begin keeping accounting information for later printing. 

Lines 9-13 create empty files (assuming they don't already exist) with the file names as shown. 

Lines 14-18 link these files with the indicated Fortran unit numbers. The files are thus the PLOT3D XYZ 
file, the PLOT3D Q file or CONTOUR plot file, the convergence history file, and the output restart flow 
field and mesh files. These lines implicitly open the files for input and output. Fortran OPEN statements 
are not used in Proteus. If the Proteus input is such that any of these files are unnecessary' (see Table 6-1), 
then the touch and In for those files can be eliminated. 

Line 19 links an existing computational coordinate system file with Fortran unit 7. If the input parameter 
NGEOM # 10, this line should be eliminated. If a restart case is being run (1REST = 2 or 3), this line 
should be replaced by the following two lines: 

In rql fort. 11 
In rxl fort. 13 

The above lines link existing input restart flow field and computational mesh files with Fortran units 1 1 and 
13. 
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Line 20 creates an executable file, p3dl0.ex, from the existing object file p3d!0.o. 

Line 21 actually runs the program. The standard Proteus input consists of all the records between line 21 
and line 22, which contains the marker 'TOD". The input could also be obtained from a file on the Cray. 
In that case, line 21 should be replaced by: 

p3dl0.ex < input 

where input is the name of the file, and line 22 should be eliminated. 

Line 23 removes the named file. 

Line 24 causes various accounting information to be printed at the end of your output. 

8.3 MODIFYING THE MASTER FILE. FULL UPDATE 


This example shows how to run with a temporarily modified version of the master file. In this particular 
case, the existing master file is modified to increase the dimensioning parameters NIP, N2P, and N3P, thus 
allowing more mesh points to be used. Since this affects almost every subroutine, the entire program is 
updated and recompiled. 


1 . 

# QSUB 

-1M 4 . OmW -IT 3600 

2. 

# QSUB 

-mb -me 

3. 

# QSUB 

-o temp. out 

4. 

# QSUB 

-eo 

5. 

# QSUB 

-r EXAMPLE 

6. 

# QSUB 

-s /bin/sh 

7. 

set -X 


8. 

ja 



9. touch plotx 

10. touch plotq 

11. touch chist 

12. touch rq2 

13. touch rx2 

14. In plotx fort. 8 

15. In plotq fort. 9 

16. In chist fort. 10 

17. In rq2 fort. 12 

18. In rx2 fort. 14 

19. In coords fort. 7 

20. cat > mods << EOF 
*id temp 

*d paramsl . 20 

parameter Cnlp = 51, n2p =51, n3p 

21. EOF 

22. update -p p3dlQ.u -i mods -c temp -f 

23. cft77 -d pq -o aggress temp.f 

24. segldr -o temp. ex temp . o 

25. temp. ex << EOD 


= 51) 


Standard Proteus input file goes here 

26. EOD 

27. rm mods temp.f temp.o temp. ex 

28. ja -cslt 


Lines 1 through 19 are the same as in the example in Section 8.2. 

Line 20 creates a Cray file called mods. The file will consist of all the records between line 20 and line 21, 
which contains the marker "EOF". Your Cray UPDATE directives and new code should therefore be in- 
serted between lines 20 and 21. The UPDATE directives and new code could also be obtained from a file 
on the Cray. In that case, lines 20 and 21 should be eliminated. 

Line 22 uses the Cray UPDATE facility to create a file called temp.f, which contains the complete Fortran 
code for the modified version of Proteus . The update command uses as input the existing Proteus program 
library p3dI0.u, and the file mods containing the UPDATE directives and new code. 

Line 23 compiles the modified program temp.f, storing the object code in the file temp.o. 
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Line 24 creates an executable file, temp. ex, from the object file lemp.o . 

Line 25 actually runs the program. The standard Proteus input consists of all the records between line 25 
and line 26, which contains the marker "EOD". The input could also be obtained from a file on the Cray. 
In that case, line 25 should be replaced by: 

temp* ex < input 

where input is the name of the file, and line 26 should be eliminated. 

Line 27 removes the named fdes. 

Line 28 causes various accounting information to be printed at the end of your output. 

8.4 MODIFYING THE MASTER FILE, PARTIAL UPDATE 

This example shows how to run with temporarily modified versions of just a few T routines. In this par- 
ticular case, the default version of subroutine INIT is replaced by a user-supplied version, and an additional 
user-supplied geometry’ option is added to subroutine GEOM. Since these changes affect only INIT and 
GEOM, only these subroutines are updated and recompiled. 

1. # QSUB -1M 4 . OmW -IT 3600 

2. # QSUB -mb -me 

3. # QSUB -o temp. out 

4. # QSUB -eo 

5. # QSUB -r EXAMPLE 

6. # QSUB -s /bin/sh 

7. set -x 

8 . ja 

9. touch plotx 

10. touch plotq 

11 . touch chist 

12. touch rq2 

13. touch rx2 

14. In plotx fort. 8 

15. In plotq fort. 9 

16. In chist fort. 10 

17. In rq2 fort. 12 

18. In rx2 fort. 14 

19. In scrl fort. 20 

20. In coords fort. 7 

21 . cat > mods << EOF 
Xid temp 
*purgedk init 
xdeck init 

U ser- supplied version of INIT goes here 

*i geom.107 

New user-supplied geometry option goes here 

22. EOF . 

23. update -p p3dlG.u -i mods -c temp -q geom init 

24. cft77 -d pq -o aggress temp.f 

25. segldr -o temp . ex temp.o p3dl0.o 

26. temp. ex << EOD 

Standard Proteus input fde goes here 

27 . EOD 

28. rm mods temp.f temp.o temp. ex 

29. ja -cslt 


Lines 1 through 20 are the same as in the example in Section 8.2. 
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Line 21 creates a Cray file called mods. The file will consist of all the records between line 21 and line 22, 
which contains the marker "EOF". Your Cray UPDATE directives and new code should therefore be in- 
serted between lines 21 and 22. The UPDATE directives and new code could also be obtained from a file 
on the Cray. In that case, lines 21 and 22 should be eliminated. 

Line 23 uses the Cray UPDATE facility to create a file called temp.f, which contains the Fortran code for 
the modified versions of subroutines GEOM and INIT. The update command uses as input the existing 
Proteus program library p3d!0.u, and the fde mods containing the UPDATE directives and new code. 

Line 24 compiles the modified versions of GEOM and INIT, contained in the file temp.f, storing the object 
code in the file temp.o. 

Line 25 creates an executable file, temp. ex, from the object file temp.o containing the modified versions of 
GEOM and INIT, and from the existing object file p3dI0.o. 

Line 26 actually runs the program. The standard Proteus input consists of all the records between line 26 
and line 27, which contains the marker "EOD". The input could also be obtained from a file on the Cray. 
In that case, line 27 should be replaced by: 

temp. ex < input 

where input is the name of the file, and line 27 should be eliminated. 

Line 28 removes the named files. 

Line 29 causes various accounting information to be printed at the end of your output. 
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9.0 TEST CASES 


In this section, two test cases are described in detail. The discussion of each test case includes a de- 
scription of the problem, listings of the Proteus input file and the JCL, and figures illustrating the computed 
results. Both cases were run using version 1.0 of the 3-D Proteus code on the NASA Lewis Cray Y-MP, 
running UNICOS 6.0, and Release 5.0.2 of CFT77. 

9.1 DEVELOPING FLOW IN A RECTANGULAR DUCT 

Problem Description 

In this test case, laminar developing flow was computed in the entrance region of a straight rectangular 
duct. The ratio of the cross-section width to height was 5:1. A schematic of part of the duct is shown in 
Figure 9.1. This flow was studied experimentally by Sparrow, Hixon, and Shavit (1967), 



Figure 9.1. Rectangular duct. 


Reference Conditions 

A convenient choice for the reference length L r was the cross-section height //, which was set equal to 
1.0 ft. The default standard sea level conditions for air of 519 °R and 0.07645 lb m /ft 3 were used for the 
reference temperature and density. The specific heat ratio y r was set to 1.4. Since the experiment was 
incompressible, the reference Mach number M r was set equal to 0.1 to minimize compressibility effects. In 
order to reach steady state within a relatively small number of time steps, the reference Reynolds number 
Re r was set equal to 60. This corresponds to Re D ~ 100, where D = 1.6667 ft. is the hydraulic diameter. 

Computational Coordinates 

For this problem a simple Cartesian computational coordinate system can be used. The physical 
( x-y-z ) and computational coordinates are thus in the same directions. Because the flow is sym- 

metric about both the y and z axes, the computational domain was just one quadrant of the cross-section. 
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Since Lr = H, the coordinate limits in the y and z directions were y min = 0, y max = 0.5, z min = 0, and 
z m cx = 2.5, respectively. The experimental data indicated that the flow should become fully developed within 
10 hydraulic diameters. The coordinate limits in the x direction were thus x min = 0 and x mc , = 16.667. 

Initial Conditions 

Constant stagnation enthalpy was assumed, so only four initial conditions were required. They were 
simply u — 1, v = w = 0, and p = 1 everywhere in the flow field. 

Boundary Conditions 

Similarly, only four boundary conditions were required at each computational boundary. At the duct 
inlet the total pressure was specified, the velocity u was extrapolated, and the gradients of the velocities v 
and w were set equal to zero. The inlet total pressure was calculated from the freestream static pressure and 
the reference Mach number using isentropic relations. At the duct exit, the static pressure was specified, 
and the velocities u , v, and w were extrapolated. The exit static pressure was found by trial and error in order 
to closely match the experimental mass flow rate and the average velocity. Since Proteus is a compressible 
code, the density will change slightly, causing the average velocity to vary with x. For the exit static pressure 
finally used, the average velocity varied from 0.973 at the inlet to 1.039 at the exit. At the walls of the duct 
no-slip conditions were used for the velocities, and the normal pressure gradient was set to zero. Symmetry 
conditions were used in both symmetry planes. These boundary conditions are summarized in the following 
table. 


Boundary 

Boundary Conditions 


u extrapolated, dvjd^ = dwjd £ = 0, pr = 1.007 

*«1 

u , v, and w extrapolated, p = 0.937 

*7 = 0 

du}br\ = 0, v = 0, dwldrj = dpjdr] = 0 

rj = 1 

u = v = w = 0, dpi dr] — 0 

c = 0 

5w/3C = 0, dv/dC = 0, w — 0, dpjdZ, = 0 

c=l 

u = v = w = 0, dpjd£ = 0 


Proteus JCL and Input File 

The Cray UNICOS job control language used for this case is listed below, including the namelist input. 
Jote that since the default for IVOUT(5) and IVOUT(6) are 30 and 40, they are et to zero in namelist IO 
to avoid printout of the static pressure and temperature. The £ indices specified for the printout correspond 
to the locations where velocity profiles were measured in the experiment. In namelist FLOW, the only 
specified reference conditions are MACHR and RER, since the desired values for the remaining ones are 
the same as the default values. ILAMV is defaulted, resulting in constant viscosity and thermal 
conductivity k. In na melist BC, the JBC values corresponding to derivative boundary conditions are posi- 
tive, specifying that two-point one-sided differences are to be used. In namelist NUM, the parameters 
IPACK and SQ for the £ and y\ directions are defaulted, which results in an evenly spaced mesh in those 
directions. The artificial viscosity is turned off because of the low Reynolds number. 

The input includes a brief explanation of each line, and should be compared with the detailed input 
description in Section 3.0. Note that it is assumed that Proteus has already been compiled, using the pro- 
cedure described in Section 8.0, with the dimensioning parameters NIP, N2P, and N3P equal to 101, 21, 
and 41, respectively. 

# QSUB -1M 8 . OmW -IT 7200 

# QSUB -mb -me 

# QSUB -o temp. out 

# QSUB -eo 

# QSUB -r ductS 

# QSUB -s /bin/sh 
set -x 

ja 

# Set up and link necessary input/output files 
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touch plotx 
touch plotq 
touch chist 
In plotx fort. 8 
In plotq fort. 9 
In chist fort. 10 

# Load and run Proteus 

segldr -o p3dl0.ex p3dl0.o 
p3dl0.ex << EOD 

Laminar developing flow in a rectangular duct 


Srstrt 

Send 

&io 

ivout=l ,2,3,32,46*0, 
iprt=l 0 0 00, 

iprtla=2,3,5,6,8,10,13,15,19,26 
iprt2= 1 , iprt3 = 2, 
iplot=2 , 

Send 
Sgmtry 
ngeom=l , 

xmin=0.0, xmax=16 . 667 , 
ymin=0.0, ymax=0.5, 
zmin=0.0, zmax=2.5, 

Send 
&f low 
ihstag=l , 

machr=0.1, rer=60., 
gamr=l . 4 , 

Send 

Sbc 

jbcl ( 1,1) =47 , jbcl ( 1 , 2) =41 , gbc 

jbcl (2 , 1 )=14, jbcl C 2 ,2) = 14 , 
jbcl (3,1) =22 , jbcl(3,2)=24, 
jbcl (4,1) =32 , jbcl ( 4 , 2) =34 , 
jbc2( 1 , 1 ) =42 , jbc2C 1 , 2) =42 , 
jbc2C2 , 1) =12 , jbc2(2,2)=ll , 
jbc2( 3 , 1 ) =21 , jbc2(3 ,2)=21 , 
jbc2C4 , 1)=32, jbc2C 4 , 2) =31 , 
jbc3C 1 , 1 )=42, jbc3C 1 , 2) =42 , 
jbc3(2, 1 )=12, jbc3(2 ,2)=11 , 
jbc3C3,l)=22, jbc3C 3 , 2) = 21 , 
jbc3(4 , 1 )=31 , jbc3(4 , 2) = 31 , 

Send 

Snum 

nl = l 01 , n2 = 21, n3 = 41, 
iav4e=0, iav2e=0, iav2i=0, 
ipackC 3) = 1 , 

sq C3,l) = l . , sq(3,2) = l.l. 

Send 
&time 
idtmod= 1 0 , 
idtau=5, cfl=10., 
ntime=1500 , 

Send 
&turb 
Send 
&ic 
u0 = l . , 

&end 

EOD 


Not a restart case . 


Print u, v, w, c p . 

Print first and last time levels only . 
40,47,63,91, £ indices for printout. 

Print at all rj indices, every other £ index. 
Write PLOT3Dj WHOLE plot files . 

Cartesian computational coordinates . 
x- coordinate limits . 
y- coordinate limits, 
z- coordinate limits. 

Constant stagnation enthalpy . 

Set Mr and Re r . 

Constant specific heats. 

( 1 , 1 ) = 1 . 007 , gbcl(l,2)=0 .937 , 

p r f p specified at £ — 0, 1. 
u extrapolated at £ = 0, 1. 

<5v/d£ = 0, extrapolated at £ = 0, 1. 

dwjd £ = 0, extrapolated at £ = 0, 1. 

dpfdrj = 0 at rj = 0, 1. 

du/drj = 0, u= 0 at rj = 0, 1. 

v = 0 at rj = 0, 1. 

dwjdrj = 0, w — 0 at rj = 0, 1 . 

dp/dC = 0 at £ = 0, 1. 

5«/3£ = 0, w = 0 at £ = 0, 1. 

<3v/d£ = 0, v = 0at£ = 0, 1. 
w = 0 at £ = 0, I. 

Use a 101 x 21 x 41 mesh. 

No artificial viscosity. 

Pack in the £ direction. 

Pack moderately near £ = 1. 

Recompute At every 10 steps. 

Spatially varying At. 

Limit of J,500 time steps. 

Laminar flow. 

Use u — p— 1, v = w = 0 oj initial conditions. 
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ja -cslt 


Computed Results 

For these calculations, the convergence criterion was that the Li norm of the residual for each equation 
be less than 10" 3 . After 1,500 time steps, the solution had not yet converged to this level, although the 
residuals wxre continuing to drop. The largest L 2 norm, which was for the x-momentum equation, had 
dropped from about 10 4 to 1 0“ 1 . The maximum and average residuals for the x-momentum equation had 
dropped from about 4 x 10 2 to 2 x 10 -3 , and 2 x 10 1 to 3 x 10~ 4 , respectively. Continuing the calculation 
past 1500 time steps resulted in further decreases in the residuals, but did not change the solution appre- 
ciably, so the calculation was stopped. The 1,500 iterations required 5,797 seconds of CPU time. 

In Figure 9.2, the computed static pressure coefficient is compared with the experimental results of 
Sparrow, Hixon, and Shavit (1967). In this figure the static pressure coefficient is defined as (#> — p)lq, 
where po is the inlet total pressure, p is the loci static pressure, and q = pi£ vg l 2 is the dynamic pressure. 
The abscissa is xj(DRe D ) x 10 2 , where x is the axial distance, D is the hydraulic diameter, and Re D is the 
Reynolds number based on the average velocity and hydraulic diameter. The computed values were gen- 
erated from the Proteus results assuming po = /w + q and p = Ua Vg — 1. The centerline value of the static 
pressure was used for p . The agreement between the computed and experimental results is good over most 
of the duct length. The largest differences are near the duct inlet and exit. This could be due to the variation 
of the density and average velocity with x, as discussed earlier. 



Figure 9.2. Computed static pressure distribution for rectangular duct flow. 


Figures 9.3 and 9.4 show the streamwise velocity profiles as functions ofy and z, respectively, at various 
x stations. The computed profiles are actually at the coordinate line in the computational grid that is 
nearest to the x station used in the experiment. No interpolation in x was done. The computed values 
were generated from the Proteus results assuming u cvg = 1. The agreement between the computational and 
experimental results is very good. 
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Streamwise Velocity, u/u avg 


Figure 9.3. Computed velocity profiles in x-y plane for rectangular duct flow. 



Streamwise Velocity, u/u avg 


Figure 9.4. Computed velocity profiles in x-z plane for rectangular duct flow. 
9.2 TURBULENT S DUCT FLOW 


Problem Description 


In this test case, turbulent flow in an S-duct was computed using first the Baldwin-Lomax algebraic 
turbulence model and then the Chien k-z turbulence model. The S-duct consisted of two 22.5° bends with 
a constant area square cross section. The geometry and experimental data were obtained from a test con- 
ducted by Taylor, Whitelaw, and Yianneskis (1982). 

Reference Conditions 

The default standard sea level conditions for air of 519 °R and 0.07645 lb m /ft 3 were used for the reference 
temperature and density. The specific heat ratio y r was set to 1.4. Since the experiment was incompressible, 
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the reference Mach number M r was set equal to 0.2 to minimize compressibility effects and, at the same 
time, achieve a reasonable convergence rate with the Proteus code. In the experiment, the Reynolds number 
based on the bulk velocity and the hydraulic diameter w’as 40,000. This value was therefore used as the 
reference Reynolds number Re r in the calculation. The reference length L r was set equal to 0.028658 ft. 
This value was computed from the definition of Re n where M r and Sutherland's law were used to compute 
u r and respectively. 

Computational Coordinates 

Figure 9.5 shows the computational grid for the S-duct, created using the GRIDGEN codes 
(Steinbrenner, Chawner, and Fouts, 1991). For clarity, the grid is shown only on three of the computa- 
tional boundaries, and the points have been thinned by a factor of two in each direction. The boundary 
grids were first created using the GRIDGEN 2D program. The 3-D volumetric grid was then generated 
from the boundary grids using GRIDGEN 3D. 26 



Figure 9,5. S-duct computational grid. 

The computational grid extended from 7.5 hydraulic diameters upstream of the start of the first bend, 
to 7.5 hydraulic diameters downstream of the end of the second bend. The grid consisted of 81 x 31 x 61 
points in the and £ directions, respectively. Since the S-duct is symmetric with respect to the rj = 1 
plane, only half of the duct was computed. To resolve the viscous layers, grid points were tightly packed 
near the solid walls using the default packing option in GRIDGEN 2D. At the gnd point nearest the w'all, 
the value of was about 0.5. 

Initial Conditions 

The computations were done in two separate major steps: a calculation using the Baldwin-Lomax tur- 
bulence model and a calculation using the Chien k-c model. To start the Baldwin-Lomax calculations, the 
default initial profiles specified in subroutine I NIT were used. Thus, the static pressure p was set equal to 
1.0, and the velocity components w, v, and w were set equal to 0.0 everywhere in the duct. To start the 
Chien k-z calculations, the initial values of u , v, w, p, and the turbulent viscosity p t w T ere obtained from the 
Baldwin-Lomax solution. The initial values of k and z were obtained using the default KEINIT subroutine 
in Proteus . 

Boundary Conditions 

For both calculations, constant stagnation enthalpy was assumed, eliminating the need for solving the 
energy equation. Therefore, only four boundary conditions were required for the mean flow at each com- 
putational boundary. In addition, for the Chien calculation, boundary conditions were required for k and 
z at each computational boundary. 


26 The form of the unformatted file created by GRIDGEN 3D was modified slightly to match that required by 
Proteus , as described in Section 3.2. 
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For the Baldwin-Lomax calculation, at the duct inlet the total pressure was specified, the gradient of u 
was set equal to zero, and the velocities v and w were set equal to zero. The inlet total pressure was cal- 
culated from the freestream static pressure and the reference Mach number using isentropic relations. At 
the duct exit, the static pressure was specified, and the gradients of u, v, and w were set equal to zero. The 
exit static pressure was found by trial and error in order to match the experimental mass flow rate. At the 
walls of the duct no-slip conditions were used for the velocities, and the normal pressure gradient was set 
to zero. Symmetry conditions were used in the symmetry plane. These boundary conditions are summa- 
rized in the following table. 


Boundary 

Boundary Conditions 

£ = o 

dujd£ = 0, v = w = 0,p T = 1.028281121 

{-1 

dufdt = dvldt = dwldt = 0, p = 0.98415512 

,7 = 0 

u = v = w = 0, dpi dr} = 0 

rj = 1 

dujdrj — 0, v = 0, dwjdrj = dpi dr] = 0 

c = o 

u = v = w = 0, dpjd£ = 0 

c=l 

u = v = w = 0, dp\d\ = 0 


For the Chien calculation, the boundary conditions for the mean flow were the same as for the 
Baldwin-Lomax calculation, with one exception. At the duct exit, the value of the static pressure was 
changed slightly, from 0.98415512 to 0.98473857, again in order to match the experimental mass flow rate. 
For the k-z equations, at the upstream boundary the gradients of the turbulent kinetic energy k and the 
turbulent dissipation rate z were set equal to zero for the first 20 time steps. After that time, the values of 
k and z were kept constant. At the downstream boundary, the gradients of k and z were set equal to zero. 
No-slip conditions were used at the surface of the plate, and symmetry conditions were used at the sym- 
metry boundary. The boundary conditions for k and z are summarized in the following table. 


Boundary 

Boundary Conditions 

o 

II 

dkld£ = dzldt = 0forn<21 
k = k 21 , z = e 21 for n > 21 

f-1 

dkldt = dtjdZ = 0 

>7 = 0 

k-z = 0 

rj = 1 

dkjdrf = dzldr, = 0 

c-o 

dkjdt; = dzjdC = 0 

C = 1 

dk/d C = dzjdC = 0 


JCL and Input for Baldwin-Lomax Calculations 

The first set of calculations used the Baldwin-Lomax algebraic turbulence model. The Cray UNICOS 
JCL file for the run is listed below. The contents of the input section of this listing should be compared 
with the detailed input description in Section 3.0. Note that for the first run, the geometry file for the S-duct 
must be available for Proteus to read in. Note also that it is assumed that Proteus has already been com- 
piled, using the procedure described in Section 8.0, with the dimensioning parameters NIP, N2P, and N3P 
equal to 81, 31, and 61, respectively. 

Since the flow field is impulsively started from zero velocity everywhere, large CFL numbers specified 
at the very beginning of the calculation might result in an unphysical flow field and cause the calculation 
to blow up. Therefore, the calculations were run with CFL = 1 for the first 100 iterations, CFL = 5 for the 
next 200 iterations, and CFL = 10 for the remaining iterations. 

# QSUB -1M lO.OmW -IT 14400 

ft QSUB -mb -me 

ft QSUB -o temp. out 

ft QSUB -eo 

ft QSUB -r sdturb 

ft QSUB -s /bin/sh 

set -x 

ja 

ft Set up and link necessary input/output files 
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touch plotx 


touch plotq 
touch chist 
touch rql 
touch rxl 

In sdcoords fort, 7 
In plotx fort. 8 

In plotq fort. 9 

In chist fort. 10 

In rql fort. 12 

In rxl fort. 14 

# Load and run Proteus 

segldr -o p3dlQ.ex p3dl0.o 
p3dl0.ex << EOD 

Turbulent $-Duct Flow, Re = 40,000, 
Srstrt 
irest=l , 

Send 

fiio 

ivout=13, 49*0, 
iprt=100000, 

iprtl=5, iprt2a=l , iprt3=2 

iplot=2 , 

Send 
Sgmtry 
ngeom=l 0 , 

Send 
Sf low 
ihstag=l , 
lr=0. 028657852, 
machr=0.2, rer=40000., 
gamr=l . 4, 

Send 

Sbc 

jbcl (1 , 1 )=47 , jbcl C 1 ,21=41 , gbclC 

jbcl ( 2 , 1 ) = 12 , jbcl ( 2,21=12 , 
jbcl (3, 1 1=21 , jbcl C3 , 2) =22 , 
jbcl(4, 11=31 , jbcl(4,2)=32, 
jbc2Cl,l)=42, jbc2( 1 , 2) =42 , 
jbc2(2, 11=11 , jbc2C 2 ,21=12, 
jbc2(3, 1 1=21 , jbc2(3 , 2) =21 , 
jbc2(4, 1 )=31 , jbc2C4 ,21=32 , 
jbc3( 1 , 1 ) =42 , jbc3( 1,2 ) =42 , 
jbc3C2,l) = ll , jbc3C2 , 2) = 11 , 
jbc3(3, 1 ) =21 , jbc3C3,2)=21, 
jbc3(4, 1 )=31 , jbc3C 4 , 2) = 31 , 

Send 

Snum 

nl=81 , n2=31, n3=61. 

Send 
Stime 
idtmod=l , 
ntseq=3 , 

idtau=5, cfl-1 . ,5. ,10 . , 
ntime=100 ,200,700, 

icteSt=4, eps=4*l . 0e~6 , 

Send 
Sturb 
iturb=l , 


Baldwin-Lomax, NASA CR 3550 

Write restart files . 


Print every 100,000 time levels. 

Print at every 5th. £ and 2nd. £ index , 
in symmetry plane only. 

Write PLOT3D plot files. 

Get computational coordinates from file. 

Constant stagnation enthalpy. 

Set L r . 

Set M r and Re r . 

Constant specific heats. 

.,11*1. 028281121, gbcl Cl ,2)=0 . 98415512, 

p Tf p specified at £ = 0 , 1 . 
dujd£ = 0 at £ = 0 , 1 . 
v = 0 f dvjd£ = 0 at £ = 0 , 1 . 
w — 0, dwjd£ = 0at £=0, 1 . 
dp/dt] = 0 at rj = 0 , I. 
u = 0 , dujdy] = 0 at rj = 0 , 1 . 
v =s 0 at ri = 0 , I. 
w = 0 , dwjdi 7 = 0 at t] = 0 , 1 . 
dp\d£ = 0 at £ = 0 , 1 . 
u=Q at £ = 0,1. 
v = 0 at £ = 0 , 1 . 
w = 0 at £ — 0 , 1 . 

Use an 81 x 31 x 61 mesh. 


Recompute Ax every time step. 

Use 3 time step sequences. 

Spatially varying Ax, CFL = 1, then 5 , then 10. 
100 time steps in first sequence, 200 in second , 
700 in third. 

Use R avi =10 ~ 6 as convergence criteria. 
Turbulent flow, Baldwin-Lomax model. 
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Send 

8ic 

Send 

EOD 

ja -cslt 

Two more runs were made with the Baldwin-Lomax model. The JCL and input for the second run 
were similar to that for the first run. The only differences were: 

1. The touch and In commands for the restart files were changed to: 

touch rq2 
touch rx2 
In rql fort. 11 
In rxl fort. 13 
In rq2 fort. 12 
In rx2 fort. 14 

2. In namelist RSTRT, IREST was set equal to 2. 

3. In namelist GMTRY, NGEOM was omitted, since the computational mesh is read from the restart 
file. 

4. In namelist TIME, NTSEQ was omitted, and the parameters NTIME and CFL were changed, as 
follows: 

idtau=5, cf 1 = 10 . , Spatially varying At, CFL= 10. 

ntime=l 500 , 1500 time steps. 

The JCL and input for the third run were similar to that for the second run. The only differences were 
that the touch and In commands for the restart files were changed to: 

touch rq3 
touch rx3 
In rq2 fort. 11 
In rx2 fort. 13 
In rq3 fort. 12 
In rx3 fort. 14 


JCL and Input for Chien Calculations 


The second set of calculations used the Chien turbulence model. It was run in a series of steps. 
The first run was a restart from the Baldwin-Lomax calculation. In this run, the gradients of k and t were 
set to zero at the inlet. Again, a small CFL number was used at the beginning of the calculation. The first 
run used CFL= 1, for just 20 iterations. The second run used CFL= 1 for the first 100 iterations, 
CFL= 5 for the next 500 iterations, and CFL =10 for the remaining iterations. Subsequent runs used 
CFL = 10 for all iterations. 


The Cray UNICOS JCL file for the first run is. listed below. The contents of the input section of this 
listing should be compared with the detailed input description in Section 3.0. 

# QSUB -1M lO.OmW -IT 300 

# QSUB -mb -me 

# QSUB -o temp. out 

# QSUB -eo 

# QSUB -r sdturb 

# QSUB -s /bin/sh 
set -x 

ja 

# Set up and link necessary input/output files 

touch plotx 
touch plotq 
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touch chist 
touch rx4 
touch rq4 
In plotx fort. 8 
In plotq fort. 9 
In chist fort. 10 
In rq3 fort. 11 
In rx3 fort. 13 
In rq4 fort. 12 
In rx4 fort. 14 

# Load and run Proteus 

segldr -o p3dl0.ex p3dl0.o 
p3dl0.ex << EOD 

Turbulent S-Duct Flow, Re * 40,000, 
Srstrt 
irest=3, 

Send 

Sio 

ivout=13, 49*0, 
iprt=100000 , 

iprtl=5, iprt2a=l, iprt3=2 

iplot=2, 

Send 
Sgmtry 
ngeom=10 , 

Send 
Sf low 
ihstag=l , 
lr=0. 028657852, 
machr=0.2, rer=40000., 
gamr=l . 4 , 

Send 

Sbc 

jbctl (1,1)= 2, jbctl ( 1,2)= 2, 
jbctl(2,l)=12, jbctl ( 2 , 2) =12, 
jbct2( 1 ,13= 1, jbct2( 1,2)= 2, 
jbct2(2,l)=ll, jbct2(2 ,2)=12, 
jbct3( 1 , 1 ) = 1, jbct3C 1 , 2) = 1, 
jbct3(2,l)=ll, jbct3C2 ,2)=11 , 
jbcl (1,1 ) =47 , jbcl C 1 , 2)=41 , gbclC 

jbcl (2, 1 ) =12, jbcl ( 2 , 2)=12 , 
jbcl C3, 1 )=21 , jbcl(3, 2)=22, 
jbcl (4 , 1 )=31 , jbclC4,2)=32, 
jbc2( 1,1) =42, jbc2( 1 , 2)=42 , 
jbc2C2,l)=ll , jbc2(2, 2)=12 , 
jbc2C 3, 1 )=21 , jbc2( 3 , 2) = 21 , 
jbc2(4,l)=31 , jbc2C 4 , 2)=32 , 
jbc3C 1 , 1 ) =42 , jbc3( 1 ,2)=42 , 
jbc3(2, 1 )=11 , jbc3C2,2)=ll , 
jbc3C3 , 1 )=21 , jbc3(3,2)=21, 
jbc3C4, 1)=31, jbc3(4,2)=31 , 

Send 

Snum 

nl=81 , n2=31 , n3=61. 

Send 
Stime 
idtmod=l , 
idtau=5 , cf 1=1 . , 
ntime=20 , 

ictest=4, eps=4*l . 0e-6 , 


Chien, NASA CR 3550 

Read restart files , previous run used Baldwin- Lomax, 


Print every 100,000 time levels . 

Print at every 5th, Z and 2nd . £ index , 
in symmetry plane only , 

Write PLOT 3D plot files. 


Get computational coordinates from file . 


Constant stagnation enthalpy. 
Set Lr. 

Set M r and Re r . 

Constant specific heats. 


dk/dZ = 0 at Z = 0, L 
dtjdZ = 0 at Z = 0, 1. 

k = 0 , dk/dri = 0 at rj = 0 , 1 . 
e = 0 , dsjdrj = 0 at 7} = 0 , 1 . 

= 0 at C = 0, 1. 
e = 0 at Z = 0, 1 . 

1,1 >=1.028281121, gbclCl,2)=0 .98473857, 

p Tr p specified at Z = 1. 

dujdZ = 0a^ = 0, 1. 
v = 0, dvjdZ = 0 atZ = 0, 1. 
w = 0, dw/dZ =0a/f = 0, 1. 

dpi dr] = 0 at rj = 0 , 1. 
u — 0, dujdrj = 0 at y\ = 0 , 1. 
v = 0 at r\ = 0, 1. 
w = 0, dwjdrj = 0 at rj = 0, 1. 
dpjdC = 0 at C = 0 , 1. 
w=0arC = 0, 1. 
v = 0 at C = 0, 1. 
w = 0 at £ — 0 , 1. 


Use an 81 x 31 x 61 mesh . 

Recompute At every time step. 

Spatially varying At. 

20 time steps 

Use R m = 10 * 6 as convergence criteria. 
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Turbulent flow , Chien model 


Send 
&turb 
iturb=20 , 

Send 
Sic 
Send 
EOD 

ja -cslt 

The JCL and input for the second run were similar to that for the first run. The only differences w T ere: 

1. The touch and In commands for the restart files were changed to: 

touch rq5 
touch rx 5 
In rq4 fort. 11 
In rx4 fort . 13 
In rq5 fort. 12 
In rx5 fort. 14 

2. In namelist RSTRT, IREST was set equal to 2, since the previous run used the two-equation tur- 
bulence model. 

3. In namelist BC, the value of JBCT1(1,1) was changed from 2 to 0, and the value of JBCT1(2,1) 
was changed from 12 to 10. This changed the boundary conditions on k and e at the upstream 
boundary from dk/dZ = d s/3£ = 0 to k = k imt and z = where k mit and z mit are the values of k and 
e at the end of the previous run. 

4. In namelist TIME, the parameter NTSEQ was added, and the parameters NTIME and CFL were 
changed, as follows: 

ntseq-3. Use 3 time step sequences . 

idtau=5, cf 1 = 1 . ,5. ,10 . , Spatially varying At, CFL= 1, then 5, then 10 . 

nti me =100,500, 400, 100 time steps in first sequence , 500 in second, 

400 in third . 

5. In the QSUB statements, the CPU time limit was raised to 14400. 

The JCL and input for the third run were similar to that for the second ru~. The only differences were: 

1 . The touch and In commands for the restart files were changed to: 

touch rq6 
touch rx6 
In rq5 fort . 1 1 
In rx5 fort . 13 
In rq6 fort. 12 
In rx6 fort . 14 

2. In namelist TIME, the parameter NTSEQ was omitted, and the parameters NTIME and CFL were 
changed, as follows: 

idtau=5, cfl=10. , Spatially varying At, CFL= 10. 

ntime=1500, 1500 time steps. 


Computed Results 


For these calculations, the convergence criterion was that the average residual for each equation be less 
than 10" 6 . However, both calculations were stopped before reaching this level of convergence when ex- 
amination of several flow-related parameters indicated that the solution was no longer changing appreciably 
with time. The average residual at the end of the Baldwin-Lomax calculation ranged from 10 -3 for the 
^-momentum equation to 3 x 10" 5 for the continuity equation. For the Chien calculation the values were 
3 x 10" 4 for the x-momentum equation and 5 x 10" 6 for the continuity equation. For both cases the resi- 
duals were continuing to drop when the calculations were stopped. A total of 4,000 iterations was used for 
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the Baldwin- Lomax calculations, requiring 31,803 seconds of CPU time. An additional 2,520 iterations 
were used for the Chein calculations, requiring 22,048 seconds of CPU time. 

In Figure 9.6, the computed flow field from the Chien calculation is shown in the form of total pressure 
contours at five stations through the duct. (The upstream and downstream straight sections are not shown.) 
As the flow enters the first bend, the boundary layer at the bottom of the duct initially thickens due to the 
locally adverse pressure gradient in that region. In an S-duct, the high pressure at the outside (bottom) of 
the first bend drives the low energy boundaiy layer toward the inside (top) of the bend, while the core flow 
responds to centrifugal effects and moves toward the outside (bottom) of the bend. The result is a pair of 
counter-rotating secondary flow vortices in the upper half of the cross-section. These secondary’ flows cause 
a significant amount of flow distortion, as shown by the total pressure contours. 

In the second bend, the direction of the cross-flow pressure gradients reverses, making the pressure 
higher in the upper half of the cross-section. However, the flow enters the second bend with a vortex pat- 
tern already established. The net effect is to tighten and concentrate the existing vortices near the top of the 
duct, in agreement with classical secondary flow theory. The resulting horseshoe-shaped distortion pattern 
at the exit of the second bend is typical of S-duct flows. 



Figure 9.6. Computed total pressure contours for turbulent S-duct flow. 

In Figure 9.7, the calculated wall pressure distribution is compared with the experimental data of Taylor, 
Whitelaw, and Yianneskis (1982). The agreement is very good. Both turbulence models correctly predicted 
the pressure trend and the pressure loss along the duct. The r and z coordinates noted in the legend are the 
same as those defined by Taylor, Whitelaw, and Yianneskis. 
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Streamwise Distance, D h 

Figure 9.7. Computed surface static pressure distribution for turbulent S-duct flow. 

In Figure 9.8, the experimental and calculated velocity profiles in the symmetry plane are shown for the 
five streamwise stations that were surveyed in the experiment. These surve r stations are at the same lo- 
cations as the total pressure contours shown in Figure 9.6. The agreement between computation and ex- 
periment is excellent for both turbulence models. The asymmetry in the velocity profiles due to the pressure 
induced secondary' motion is correctly predicted. 
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Figure 9.8. Computed streamwise velocity profiles for turbulent S-duct flow. 
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